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The subject of cosmological hydrogen recombination has received much attention recently because 
of its importance to predictions for and cosmological constraints from CMB observations. While 
the central role of the two-photon decay 2s — > Is has been recognized for many decades, high- 
precision calculations require us to consider two-photon decays from the higher states ns,nd — > Is 
(n > 3). Simple attempts to include these processes in recombination calculations with an effective 
two-photon decay coefficient analogous to the 2s decay coefficient A2 S = 8.22 s _1 have suffered from 
physical problems associated with the existence of kinematically allowed sequences of one-photon 
decays, e.g. 3d — > 2p — > Is, that technically also produce two photons. These correspond to reso- 
nances in the two-photon spectrum that are optically thick to two-photon absorption, necessitating 
a radiative transfer calculation. We derive the appropriate equations, develop a numerical code to 
solve them, and verify the results by finding agreement with analytic approximations to the radia- 
tive transfer equation. The related processes of Raman scattering and two-photon recombination 
are included using similar machinery. Our results show that early in recombination the two-photon 
decays act to speed up recombination, reducing the free electron abundance by 1.3% relative to the 
standard calculation at z = 1300. However we find that some photons between Lya and Ly/3 are 
produced, mainly by 3d — > Is two-photon decay and 2s — > Is Raman scattering. At later times 
these photons redshift down to Lya, excite hydrogen atoms, and act to slow recombination. Thus 
the free electron abundance is increased by 1.3% relative to the standard calculation at z = 900. 
Our calculation involves a very different physical argument than the recent studies of Wong & Scott 
and Chluba & Sunyaev, and produces a much larger effect on the ionization history. The implied 
correction to the CMB power spectrum is neligible for the recently released WMAP and ACBAR 
data, but at Fisher matrix level will be 7a for Planck. 

PACS numbers: 98.70.Vc, 98.62.Ra, 32.80.Rm 



I. INTRODUCTION 

The anisotropics of the cosmic microwave background 
(CMB) are one of the most important tools for cosmol- 
ogy. The power spectrum from the Wilkinson Microwave 
Anisotropy Probe (WMAP) satellite has enabled preci- 
sion measurement of the composition of the high-redshift 
universe (the baryon/photon and dark matter/baryon ra- 
tios), the distance to the surface of last scattering, and 
the primordial power spectrum (including its spectral in- 
dex n s ) [H, • The first two of these have proven essen- 
tial to studies of dark energy, e.g. by breaking the ap- 
proximate fl m — w degeneracy in low-redshift supernova 
data, while the mesurements of n s have begun to rule 
out interesting classes of inflationary models. This trend 
can be expected to continue with the upcoming Planck 
satellite and several ground- and balloon-based experi- 
ments to measure small-scale temperature anisotropics 
and CMB polarization. In particular, Planck should be 
able to map the entire CMB sky with signal-to-noise ra- 
tio greater than 1 down to scales of I ~ 1600, which in 
principle enables measurement of n s with an uncertainty 
of - 0.003 ||. 

While the raw sensitivity of the future CMB projects is 
spectacular, there are significant challenges for both ex- 
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pcrimcnt and theory. It is well-known that these projects 
require removal of instrumental signatures, foregrounds, 
and secondary anisotropics, and highly accurate beam 
maps. In contrast, the theory of primary anisotropics 
is generally regarded as simple: it is linear perturbation 
theory whose development was completed more than a 
decade ago, and which can now be solved by public com- 
puter codes that agree to within 0.1% 0]. However, the 
linear perturbation theory calculation requires knowledge 
of the free electron density (and hence Thomson opacity) 
in the unpertubed Universe. This free electron density is 
determined by the complicated non-equilibrium physics 
of recombination, as first noted forty years ago by Pee- 
bles [H and Zel'dovich et al. @ . For this reason, a great 
deal of effort in the 1990s was aimed at solving cosmic 
recombination. This culminated in the recfast code by 
Seager et al. which with some recent improvements 

@ is used to compute the recombination history in most 
of today's CMB prediction codes. The Seager et al. anal- 
ysis is based on the "multi-level atom" (MLA) method, 
in which one writes down the occupation probabilities for 
each level of each atom or ion, and then constructs a set 
of evolution equations that includes radiative transition 
rates, photoionization and radiative recombination rates, 
and collision terms. A separate allowance is made for the 
two-photon transition H I 2s — > Is, and one must also in- 
clude an equation for the matter temperature (including 
all heating and cooling terms) since at late times this is 
not in equilibrium with the CMB temperature. The H I 
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and He II Lyman-series lines, and the He I n}P° — l 1 5o 
lines, become optically thick and are treated using the 
Sobolev [l(| approximation. 

In the past several years, there has been a resur- 
gence of interest in the recombination problem, driven 
primarily by the new CMB experiments [1 EH and 
also the possibility of measuring the spectral distortion 
El El El El EE El Gl- Several groups have added 
additional physics to the Seager et al. MLA. This in- 
cludes two-photon transitions from higher excited states 
of H I (us, nd -> Is) El H HI and He I H3; stimulated 
two-photon decays and two-photon absorption [H [24[ ; 
separate consideration of each nl sublevel of H I instead 
of assuming a statistical distribution of I 12511; semifor- 
bidden and forbidden transitions in He 1 [20l l26l , [27j ; 
non-Sobolev behavior in lines with significant continuum 
opacity and partial redistribution (2a . |28| ; and recoil of 
H 1 atoms in Lya scattering p9| . 

This paper is concerned with the two-photon transi- 
tions during the H 11— >H 1 recombination, with special 
attention given to the ns, nd — > Is decays with n > 3. 
The physics of the two-photon process has been well- 
understood since the work of Goppert-Mayer [3(j, and 
successively more accurate computations of the 2 s life- 
time have been made in the following years iil. 
The importance of the 2s — ► Is decay for cosmological 
recombination was recognized in the 1960s [HE and it is 
typically included in MLA codes by adding a 2s — > Is de- 
cay with a decay rate A2 S = 8.22 s _1 and a thermal two- 
photon absorption rate set by the principle of detailed 
balance. However, cosmologists paid little attention to 
the two-photon decays from other states until the recent 
work of Dubrovich & Grachev E^l (hereafter DG05). De- 
spite the major recent advances in this subject [2fl. [2l|. 
we still lack a fully self-consistent treatment of the prob- 
lem. 

DG05 attempted to extend the MLA to two-photon 
transitions from highly excited (n > 3) levels by defin- 
ing analogous two-photon rates A n s,nd- Unfortunately 
they discovered that this does not work: the matrix el- 
ement has poles corresponding to the "1+1" photon de- 
cays of the form ns,nd — > Np — > Is with 1 < N < n. 
This results in a very fast decay rate for the n > 3 
levels, with most of the decays producing a photon in 
a Lyman-series resonance. However since the photons 
emitted in resonance lines have already been counted in 
the one-photon treatment, DG05 recognized that to avoid 
double-counting, they should somehow remove the reso- 
nance from their two-photon rates. They did this by 
keeping only one pole in the matrix clement, that corre- 
sponding to the np intermediate state. This led them to 
very fast decay rates that scale as h. ns , n d °c n for high n 
and a several percent acceleration of hydrogen recombi- 
nation. Wong & Scott (2(j (hereafter WS07) compared 
the DG05 result for n = 3 to calcuations [H HH that 
included all non-resonant poles; they found a lower rate 
than DG05 and re-scaled DG05's A ns . Il( j values appropri- 
ately. This resulted in a maximum change of 0.4% in the 



electron abundance. Chluba & Sunyaev [2l| (hereafter 
CS08) argued that one should cut off the two-photon de- 
cay rate not by selecting poles but based on the physical 
criterion of whether one of the photons is immediately 
re-absorbed in a Lyman-series line. They computed full 
two-photon spectra and implemented a cutoff based on 
their calculation of absorption in the red damping wing of 
Lya. This procedure has the distinct advantage of being 
convergent as n max — > 00. The results do depend some- 
what on the nature of the cutoff (xr, c in CS08 notation) 
but generally give changes in the electron abundance of 
several tenths of a percent. 

In this paper, we write down the full system of equa- 
tions for two-photon transitions, including emission, ab- 
sorption, and Raman scattering. We then take two ap- 
proaches to solving these equations. We first present a 
numerical approach in which the two-photon continuum 
is discretized and turned into an effective MLA with vir- 
tual levels. We then present an analytic approach in 
which the two-photon transitions are approximated as 
effective corrections to the Lya and Ly/3 decay rates and 
added into the standard MLA with no virtual levels. The 
first approach, like most fully numerical approaches, has 
the advantage of solving the two-photon equations with 
accuracy limited only by step sizes (in both frequency and 
time!) and machine precision. Unfortunately, the results 
simply come out of the computer and are difficult to in- 
terpret. The analytic approach allows one to understand 
on paper the cancellations that cause the two-photon ef- 
fect to be finite and independent of cutoffs. It also serves 
as a check on the fully numerical result. In the future 
it may also have some value in "fast" recombination cal- 
culations that could be embedded in Markov chains for 
cosmological parameter estimation. 

This paper is not intended to be a complete solution 
of the hydrogen recombination problem. A number of ef- 
fects that are not considered here such as the very high-n 
levels, Lya diffusion and recoil, and feedback from the 
helium spectral distortion EE El, SI HI IM HI HI, 
may also be important if sub-percent accuracy is desired. 
Rather, the purpose of this paper is to clear up the phys- 
ical picture of the two-photon transitions and estimate 
their effect on recombination, as a first step toward a 
treatment that will ultimately include these other pro- 
cesses as well. 

This paper is organized as follows: in Section |TT] we 
review the multi-level atom method that is traditionally 
used in recombination studies and the sometimes-used 
steady-state approximation. In Section IIIII we describe 
how to graft two-photon transitions onto this framework, 
the problems that have occurred in past efforts to define 
a two-photon decay rate from highly excited states, and 
the resolution of these problems. We describe the nu- 
merical method for solving the problem and its results 
in Section IIV1 and an analytic approximation in Sec- 
tion |Vl The effect on the CMB power spectrum is as- 
sessed in Section IVII We conclude in Section IVIIl The 
appendices cover several technical details: Appendix [A] 



3 



considers the early-time matter temperature evolution, 
and Appendix[B]computes the two-photon decay, Raman 
scattering, and two-photon recombination rates. 

Throughout this paper, we use a background cosmol- 
ogy with n m h 2 = 0.13, Q b h 2 = 0.022, T C mb = 2.728 K, 
helium mass fraction Y = 0.24, and an effective number 
of massless neutrinos N v = 3.04. This implies the derived 
parameters /h c = 0.0795 (He:H ratio by number) and ra- 
diation density £l r h 2 = 4.196 x 10 -5 . (These parameters 
are no longer up-to-date but should suffice to illustrate 
the physics and the magnitude of the corrections due to 
two-photon transitions.) The cosmological constant is 
neglected because it has no influence on the recombina- 
tion era, £Ia(z) < 10~ 8 . Reduced matrix elements are 
normalized following the convention of Berestetskii et al. 
[39j . We will define a two- photon decay from a highly ex- 
cited (n > 3) level to be "sub- Lya" if both emitted pho- 
tons are below the Lya frequency and "super- Lya" if one 
photon is emitted above the Lya frequency, in order to 
avoid repeating these cumbersome descriptions through- 
out the paper. 



II. THE STANDARD MULTI-LEVEL ATOM 

In this section, we revisit the basic equations of the 
multi-level atom (Section III A[) and the steady-state ap- 
proximation (Section III Bp . 

Most of the notation of this section is conventional. 
We define nn to be the total physical density of hydrogen 
nuclei (units of cm -3 ) and define the abundance of hy- 
drogen atoms in level i relative to this total: x% — rii/riji- 
A similar definition is used for free electrons and protons, 
x e = n e /nn and x p = n p /nH- We denote the energy of a 
bound state by Ei < and the energy of a free electron 
by E e > 0. The degeneracies of states are denoted by 
gf, in the case of hydrogen where we do not resolve fine 
structure levels, g n i = 2(21 + 1). Einstein coefficients are 
denoted Aji for a transition from upper level j to lower 
level i, and transition energies/frequencies are denoted 



E 



Ej — Ei and Vji 



Eji/h. The photon phase space 
density is f(E). When considering the phase space den- 
sity just above (blueward of) or below (rcdward of) a 
line, we use the notation 



f Jl± =f(E Jl ±e) 



(1) 



where e is the line width, which is assumed to be infinites- 
imal (i.e. the time required to redshift through the line 
is assumed to be much less than the duration of recom- 
bination). Temperatures are quoted in energy units, so 
that Boltzmann's constant is equal to 1. 



A. Basic equations 



The evolution equation for xi is 

Xi = Xi\bb + %i\b{ + %i\2 7 , 



where we have separated the net production of level i 
into terms coming from radiative bound-bound, radiative 
bound-free, and two-photon processes. 
The bound-bound term is 



Xi\hh 



— 2_j AjiPji 



(! + fji+)xj - —fji+Xi 

9i 



+ 2_j Aij Pij 



fij + Xj (1 ~t~ fij+)Xi 

9j 



,(3) 



where the sums are over levels j that are above (j > i) 
or below (j < i) the energy of level i; Aji is the Einstein 
coefficient; gi and gj are the level degeneracies; and fji + 
is the photon phase space density on the blue (incoming) 
side of the line with frequency Vji = (Ej —Ei)/h, which is 
necessary for following absorption and stimulated emis- 
sion. Here Pji is the Sobolcv escape probability for the 
ji line, which is the probability that a photon emitted 
in the line will escape via redshifting without being re- 
absorbed. It is given by 



P, 



with the optical depth 



1 - e- T J< 



_ c"n H g 3 



SttFzA \g 



(4) 



(5) 



see e.g. Ref. [8| for a derivation. In practice the optical 
depth is Tji <C 1 and Pji ps 1 except for the Lyman-series 
lines. The jump condition for the radiation field across 
the line is obtained by finding the net rate of radiative 
de-excitations (i.e. photons emitted in the line), and mul- 
tiplying by the conversion factor 8irHvji/ (c 3 nn), which 
is the number of photon modes that redshift through the 
line per hydrogen nucleus per unit time: 



fji- — fji+ — — ^jiPj 



ji 



Xj(l + fji+) - —•'•-/.;<• 

9i 



Xi\U 



. (6) 

The bound-free term includes spontaneous and stimu- 
lated recombination, and the inverse process, photoion- 
ization (see e.g. Sec. 2.3.1 of Ref. \^). It is 

J ai (E e ){p M (E e )n H x e x p [l + f(E e - Ei)] 

where E e is the energy of the free electron (more ac- 
curately, the center-of-mass energy of the electron and 
proton), f(E e — Ei) is the photon phase space density 
at specified energy, a>i(E e ) is the recombination coeffi- 
cient to level i for an electron and proton, and Pu_(E e ) is 
the Maxwcllian probability distribution for the electron 
energy at matter temperature T m : 



(2) 



2 E 1/2 

P M (E e ) dE e = -= -5_ e -^/ T » dE P 
V 71 " T, 



3/2 



(8) 
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The 1 + f{E e — Ei) factor accounts for stimulated recom- 
bination. The factor dg/(dE e dV), i.e. the number of 
available electron states per unit volume per unit energy, 
is required in order to satisfy the principle of detailed 
balance. It is equal to 



<kl 8\Z2?r 3/2 1/2 



dE P dV 



h 3 



(9) 



where fi is the proton-electron reduced mass, \x 
m e /(l + m e /m p ). Thus we may write: 



Xi\bf 



9 f E 1 ^ 2 
x^n H x e x p [l + f(E e - Ei)]( 



— Ee/Tfn 



2(27rT m ) 



3/2 



^ E lj2^ f{Ee _ Ei) -\ dEem 
Pi ' 



h 3 - gi 

Finally the usual two-photon term is 

i2s|2 7 = -*i s |2 7 = A 2s [~x 2s + x ls e~ E2s ' ls/Tl ^ , (11) 

where A2 S is the two-photon decay rate from the 2s level, 
and the second term accounts for detailed balance. 

Note that the bound-bound radiative transition rate 
requires knowledge of the photon phase space density 
f(Eji + ) on the blue side of the line. Seager et al. [8[ 
assumed that this is simply a blackbody function, 



f{E ji+ ) 



1 



i /T m _ l - 



(12) 



In reality the phase space density for the Lyman-series 
lines will be greater than the blackbody value, because 
photons from the Ly(n + 1) line will redshift downward 
and add to the Lyn line [2y, |4(| ■ In the absence of any 
processes that operate between the Lyn and Ly(n + 1) 
lines, one may simply use phase space conservation along 
a photon trajectory: 



fnp,ls+(z) — /(n+l)p,ls- 



l-(n + l)- 
1 -n" 2 



-(! + *)-! 



(13) 

This equation refers to a photon phase space density 
at some earlier time. We generate a lookup table of 
/(n+i)p,is- as wc proceed through recombination, and 
when a value is required we obtain it by 3-point quadratic 
interpolation. This method requires that the step size 
Aa/a be smaller than the spacing of the Lyman lines, 
v (n+i)p,isl v np,is ~ 1) which imposes a significant but 
tractable computational burden. 
The matter evolution equation is 



-2HT„ 



8x e <jTa r T^ 



3(1 + /ho + x e )m e c 



(T r - T m ), (14) 



where the second term comes from electron scattering 
m. It contains the Thomson cross section ctt and the 



radiation constant a r . The second term is very effective 
at driving the matter temperature toward the radiation 
temperature at early times. Later when x e and T r 4 de- 
cline the matter temperature falls below the radiation 
temperature. Since the radiation temperature scales as 
T r oc a -1 , we may write this equation as 



dt V T r 



T r 3(1 + / Hc + x e )m e c 



B. Steady-state approximation 



(15) 



We now introduce the steady-state approximation, 
which is used to accelerate solution of the hydrogen re- 
combination problem by eliminating the need for a stiff 
differential equation solver. 

The above equations for the excited hydrogen levels 
(i =/= Is) can be written in the form 

Xi = Si + ^2 {Rji'Xj - RijXi) - PiXi - jiXi, (16) 

where Si is a source term that depends on excitations 
from Is — > i and recombinations to level i, and Rij is 
the transition rate from level i to level j, fli is the rate 
of ionizations from level i, and ji is the rate of decays 
(spontaneous, stimulated emission) from level i to the 
ground state. The coefficients {si,i?y , /3i,7i} in general 
depend on the ionization fraction, radiation and matter 
temperatures, free electron density, and spectral distor- 
tions (since they involve fji+)- The explicit expressions 
arc: for the source term. 



— — A p f 

°i ^ ^i,ls-^ i,ls Ji 7 \s-\-X\s 



+ -^=nB_x e x p J [I + f{E e - Ei)} 
x exi{E e )El' 2 e - Ee/TmdEp 



1 m 

for the transition term, 

for the ionization term, 
2 7 / 2 V /2 



(17) 



(18) 



h 3 



[ ai(E e )—f{E - Ei)dE e ; (19) 
J 9i 



and for the decay term, 

7l = Ai,uPi,uO- + fi,u+) + k2sA 2s - (20) 

In these equations, there is a dependence on x e , but 
almost no dependence on the XiS. The exceptions are 
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(i) formally x\ s = 1 — x e — Xi depends on the popu- 
lations of the excited states, and (ii) the Sobolev prob- 
abilities Pij contain a small dependence on the excited 
states. However past analyses have shown that the pop- 
ulations of the excited states are very small (the maxi- 
mum is X2 P = 5.5 x 10~ 14 at z = 1360) so we may take 
X\ a ~ 1 — x e . Also the Sobolev escape probabilities for 
transition between two excited state are always very close 
to 1 (the peak optical depth for such lines is 5.5 x 10~ 4 
for Ha at z = 1400, corresponding to an escape prob- 
ability of 0.99972). Finally the Sobolev probability for 
transitions involving the Is level are very small, but the 
high optical depth is due entirely to the Is level - the 
contribution from the excited level np changes the opti- 
cal depth by a factor of x np /3xi s which has a maximum 
value of 1.7 x 10~ 12 (for n = 2, z = 1600; at later times 
or for higher n it is lower). This leads to a convenient 
re-packaging of Eq. (fT6|) : 

X{ Si ^ ^ ^ij ! ( ^ 1 ) 

3 

where 

Tij = Sij (pi + 7i + ^ fc J _ /, '« ( 22 ) 

is a square matrix. 

The excited states of a hydrogen atom typically have 
very short lifetimes compared to the recombination 
timcscale: the longest intrinsic lifetime is that of 2s 
(A^ 1 = 0.12s), and even taking into account Sobolev 
suppression the intrinsic lifetimes of the p-states are 
short, e.g. for 2p we have (^4.2p,is-P2p,is) ~ 1 s when the 
Lya optical depth reaches its peak of ~ 6 x 10 8 . For com- 
parison, hydrogen recombination takes severalx 10 12 s. 
Therefore we expect that to very high accuracy the ex- 
cited level populations should be in steady state given 
Eq. (l2Tj) . i.e. wc should have 

Xi « ^(T- 1 )^- (23) 

3 

and 

X e ~ —Xis = — (liXi — ^Aj isPi isfi is+Xi^j . 

i>ls 

(24) 

Formally it is the minimum eigenvalue t r 1 of the matrix 
T that indicates how rapidly a steady state solution is 
approached. During the redshift range covered by our 
code the relaxation time t r (reciprocal of the minimum 
eigenvalue of T) has a maximum of 0.8 s, or ~ 10~ 12 of 
the duration of recombination. Thus the steady state 
approximation should very accurately describe the pop- 
ulations of the excited states and the ionization history 
as a function of time. (It may not accurately describe 
the low-frequency spectral distortion since the latter de- 
pends on slight deviations of the excited state ratios from 



thermal equilibrium, i.e. one is taking differences of ex- 
cited state populations that are nearly zero. In general 
the spectral distortion is much more sensitive than the 
recombination history to numerical errors [Hi EH-) 

In the case of the Lyman lines where there is significant 
feedback, we need the phase space density f( n +i)p,is-i 
which is easily obtained from Eq. (J5]) once we have ob- 
tained X( n +i) p from Eq. (|23p . 

In order to eliminate the need for a stiff ODE solver we 
also need to remove the rapidly decaying mode from the 
matter temperature equation, Eq. (|15p . At early times 
this can be done using the first-order asymptotic solution, 

_ 1 _ 3(1 + / He + x e )m e cH 

(See Appendix |A"1 for a derivation). This equation can be 
used to solve algebraically for T m at each step in the inte- 
gration, rather than requiring a stiff integrator. In fact it 
is quite accurate at early times: the recfast code, which 
fully integrates Eq. (fT5|) at low redshifts, gives a value for 
T m differing by < 0.01% from Eq. ^ for all z > 700. 
(At early times, recfast sets T m = T r since the Comp- 
ton equilibrium time is fast, but we have checked that at 
z < 850, recfast has switched on the integration of the 
T m equation so this is a fair comparison.) Accurate calcu- 
lations at lower redshift would require one to "switch on" 
the differential equation for T m once the matter-radiation 
equilibration time (HJ)^ 1 becomes a significant fraction 
of the recombination time; at this point it would be pos- 
sible to follow T m with a non-stiff integrator. 

In most implementations of the ML A, one first groups 
states together and then solves either the "full" Eq. (|2|) 
or the steady state problem, Eq. For example the 

original calculation of Peebles [5j used a single, lumped 
excited state and the steady state approximation (with 
some additional approximations such as neglect of stim- 
ulated recombination). 

III. EXTENSION TO TWO-PHOTON 
TRANSITIONS 

Now we would like to extend our analysis to include 
two-photon decays from the high- lying levels. The selec- 
tion rules for two-photon electric dipole decays restrict 
us to excited levels ns and nd (two-photon np — > Is de- 
cays are parity-forbidden). At first sight this is a simple 
problem. The two-photon decay rates can be calculated 
in tree-level QED, as done in Appendix 151 and the total 
transition rate A„; can be found using Eq. (|B6|) . This 
rate can be added to the evolution equations in analogy 
to Eq. (nj. 

Unfortunately this program does not work because 
the desired rates of e.g. spontaneous two-photon decay 
3d — > Is double-count the "1+1" or cascade decay pro- 
cess 3d — > 2p — > Is, and do not correct for re-absorption 
of the emitted photons. Wc develop our understanding of 
this issue by first discussing the two-photon decay rates 
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in Sec. IIII Al and then reviewing previous attempts to in- 
corporate two-photon transitions in the MLA f Sec. IIII Bj) . 
We discuss the physical resolution of the double-counting 
and re-absorption problems in Sec. IIII CI We then write 
down the full radiative transfer equations in Sec. IIII Dl 
The principal result of this paper will be the solution of 
this last set of equations. 

A. Rate coefficients 

The theory and computation of two-photon decay rate 
coefficients is reviewed in Appendix [B] Here we merely 
recount the most important facts and results. The two- 
photon decay is 

H(raZ) -> H(ls) +hv+ hv', (26) 

where we define v to be the higher-energy photon and v' 
to be the lower-energy photon, so that there is no issue of 
indistinguishable particles and the corresponding factors 
of 1/(2!). The two photons have frequencies that sum to 
(1 — n~ 2 )7Z due to conservation of energy The rate of 
decays is given by an integral of the differential rate 

r„^i s = / — dv, (27) 

J(l-n- 2 )Tl/2 UV 

where in a thermal radiation bath the actual decay rate 
is related to the spontaneous rate dA/dv by 

fE = ^(1- e- h "/ T ')- 1 (l - e-^'^)- 1 . (28) 
dv dv 

The increase over the vacuum (T r = 0) rate is due to 
stimulated two-photon decays [23j |. 

Panel (a) of Fig. [1] shows the two-photon decay rate 
from 2s — * Is at several temperatures. The differential 
decay rate dT / dv is well-defined and integrable at all tem- 
peratures, and this decay presents no special problems. 

Panel (b) of Fig. [1] shows the 3s — * Is decays. Here 
Eq. (|2"8"|) contains a resonance at v = 37Z/4, i.e. the 
Lya frequency. This corresponds to the pole in the ma- 
trix element from the 2p intermediate state. Taking the 
two-photon decay formula at face value the total rate 
r 3s ^i s is infinite due to the ~ 1/Av 2 divergence in the 
rate. If one includes the imagniary energy ("pole dis- 
placement") of 2p due to its finite lifetime (1.6 ns) then 
one gets a finite decay rate, but very large: in vac- 
uum r3 S ^i s = 6.3 x 10 6 s _1 . This is equal to the one- 
photon transition rate A 3s ^p (see discussion in CS08) for 
a simple physical reason: the sequential or "1+1" decay 
3s — > 2p — > Is is simply a two-photon decay with the 
photons on resonance. (The one- and two-photon con- 
tributions can be separated as different diagrams that 
contribute to the width or imaginary energy shift of the 
3s level SsAE 3s [UEi], but 3A£ 3s does not give us the 
frequency distribution of emitted photons, which we need 
in our calculations.) A similar problem occurs for the 3d 
initial level, as seen in panel (c) of Fig. [T] 



A similar problem occurs for Raman scattering: 

H(nZ) + hv' -> H(ls) + hv. (29) 

This time it is the difference of frequencies that is fixed 
by energy conservation: v = v' + (1 — n~ 2 )lZ. There is 
no Raman scattering in vacuum (T r = 0) as it requires 
an initial-state photon, but one can define a Raman- 
scattering rate in analogy to Eq. (|28p at finite tempera- 
ture. This is shown for the 2s initial state in the panel 
(d) of Fig. [TJ In this case there are resonances when the 
virtual energy of the intermediate state corresponds with 
an energy level of hydrogen: at v = 81Z/9, 157^/16, etc., 
and they reproduce the thermal rate for the sequential 
1+1 photon transitions 2s — * 3p — > Is, 2s —tAp—t Is, 
etc. 

For completeness, we also discuss two-photon recom- 
bination, 

H+ + e~ ->R(ls) + hv + hv'. (30) 

The differential cross sections da/dv for this process are 
shown in Fig.[5J It too exhibits Lyman-series resonances. 
Note that the differential rate coefficient a^(v, v') is 
equal to the incident electron velocity times da/dv: 
a^ 2 '(v, v') = v e da/dv. 

At this stage, we have two basic problems. One is 
the double-counting problem: we have already included 
1+1 photon transitions in the MLA, so it seems incor- 
rect to recount the resonant two-photon transitions. The 
other is the re-absorption problem: the very large two- 
photon transition rates near resonance imply large op- 
tial depth for re-absorption of some of the emitted radi- 
ation, combined with the possibility that the Lyman-line 
spectral distortions could be immediately re-absorbed. 
The next section critically discusses previous attempts 
to solve these problems. 

B. Past attempts at inclusion in MLA codes 

The double-counting and re-absorption problems were 
recognized even in the very first attempts to include two- 
photon decays in the MLA. There have been three major 
approaches suggested in the recombination literature so 
far: those of DG05, WS07, and CS08. Here we describe 
each of them and explain why, while each of these meth- 
ods contains important physical insights, none of them 
are both complete and applicable to the problem at hand. 
In each of these cases the authors themselves noted that 
their solution was only a first approximation and that 
ultimately a radiative transfer calculation was required; 
in particular CS08 provided a detailed description of the 
problem. 

DG05 avoided the double-counting problem by remov- 
ing the 1+1 poles from the two-photon decay matrix ele- 
ment M.. For example, in the case of the id — > Is decay 
they removed the 2p intermediate state. They then ar- 
gued that in ns, nd — > Is decays, the largest contribution 
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(a) 2s-1s two-photon decay rate 



(b) 3s-1s two-photon decay rate 
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(d) 2s-1s Raman scattering rate 
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FIG. 1: The two-photon transition rates for the most important excited levels of H I. The horizontal axis shows the frequency 
of the higher-energy photon (for 2y decays) or the outgoing photon (for 2s — > Is Raman scattering). The differential transition 
rate dY jdv is shown at several different temperatures and in vacuum (T r = 0). The two-photon spectra are shown from half the 
maximum frequency (i>Wx/2) to ^ max and are symmetric. The Lya resonance at v = 0.757?. is easily seen in the 3s, 3d — ► Is 
rates. 
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FIG. 2: The differential cross section da jdv for two-photon 
recombination at three incident electron energies (0.1, 0.05, 
and 0.025M?). 



to the matrix clement comes from the np intermediate 
state. Keeping only this term in M. they arrived at a 



a DG05 
Ly nl 



21 + 1 



| (nZ||r| \np) (np\ |r| |ls) 



(31) 



which they normalize to the 2s recay rate A2 S = 8.22 s _1 . 
They then substitute this rate coefficient into their MLA 
code in the same way that the traditional treatment in- 
serts A-2s- The DG05 rate coefficients scale as oc n and 
so a cutoff must be imposed; they impose such a cutoff 
when the wavelength of emitted radiation is less than the 
size of the atom (~ a$n 2 ) because at higher n the electric 
dipolc formula is no longer valid. 

This approach obviously eliminates the double- 
counting problem. Also by getting rid of the large dif- 
ferential rate dA/dv near resonance it eliminates the re- 
absorption problem. However it ignores the destructive 
interference in M. that arises from consideration of the 
different intermediate states; it predicts dAjdv oc n for 
large n and fixed v, whereas the actual scaling when in- 
terference is considered is oc n~ 3 . This occurs because 
while the largest contribution to the ns,nd — > Is ma- 
trix element conies from the np intermediate state, the 
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neighboring intermediate states (n ± l)p, (n ± 2)p, etc. 
give opposite sign contributions and result in a cancella- 
tion of the matrix element that becomes more exact as 
n — > oo [22j| . Thus to have even a qualitatively correct 
result for dA/dv we must include the troublesome inter- 
mediate states such as (n — l)p, even if they are far from 
resonance. 

WS07 presented an improved treatment of the prob- 
lem. They used rate coefficients for 3s — > Is and 3d — > Is 
two-photon decays that included all the n > 3 interme- 
diate states (i.e. those that cannot be on-shell) [H|, [H[ • 
From the n = 3 levels this gives a two-photon decay rate 
of 0.0664 times the DG05 estimate, a result of the 4p, 
5p, etc. intermediate states partially cancelling the 3p 
contribution to M.. WS07 did not have calculated rate 
coefficients for n > 3 so they re-scaled all of the DG05 es- 
timates by a factor of 0.0664 and included them in their 
MLA code. The WS07 method is an improvement be- 
cause it partially takes into account the destructive in- 
terference in the matrix element. Nevertheless, it still 
does not take into account all intermediate states and 
retains the incorrect dA/dv cx n scaling. 

CS08 presented the full calculation of the two-photon 
decay rate dA/dv including all intermediate states. How- 
ever by including these intermediate states they re- 
introduce the double-counting and re-absorption prob- 
lems. CS08 argued that the double-counting problem 
should be solved by defining a "1+1 decay profile" (</> ra + s 7 d 
in their notation) consisting of the sum of Lorcntzians at 
each resonance. The decay rate from the 1+1 profile 
is / <fil£s? d dv = Y^=2^nl,n' P - They then subtract this 
from the full two-photon decay profile to give a "pure" 
two-photon profile. Their effective two-photon rate is ob- 
tained by integrating this pure two-photon profile, 



A CS08 



(l-n- 2 )K/2 



dA nl 
dv 



O ) *>, 



(32) 



where v c is a cutoff frequency. (CS08 denote this rate 
by A 2 " 1 . .) The "obvious" choice for the cutoff frc- 

quency would be the maximum frequency (1 — n~ 2 )lZ, 
but CS08 argued that a lower cutoff should be used to 
solve the re-absorption problem. All photons emitted 
blueward of Lya will be re-absorbed later in recombina- 
tion. Furthermore, photons emitted rcdward of Lya will 
be re-absorbed if they are within the optically thick part 
of the red Lya damping wing, so v c should be somewhat 
below fLyc*. CS08 parameterized this in terms of the 
number of Doppler widths xr. c - 



2T r 

TOHC 2 



(33) 



(Note the sign convention that xt, c < corresponds to 
the red side of Lya.) The pure two-photon rate A^f 08 
depends only logarithmically on xr, c because the 1+1 
profile removes the order (y — v^ ya )~ 2 term from the se- 
ries expansion of the integrand in Eq. (|3"2"|) . so that at 



small v — v-^ya the leading term is of order (y — J^Lya) -1 , 
whose integral is logarithmically divergent. CS08 consid- 
ered several values of xr, c ranging from — 10 4 to — 10 5 , 
and found that their correction to the free electron abun- 
dance was < 0.5% in each case. They also found that 
while the full two-photon profile dA n i/dv must be non- 
negative, there is no restriction on its sign after the 1+1 
profile has been subtracted; their rate coefficient A^ sos is 
positive for the d initial states and negative for s initial 
states. 

While the CS08 treatment is the most complete thus 
far, there are three areas in which it could be improved. 
One is that it contains a free parameter xr, c that affects 
the results and it is not obvious which value should be 
used. A second issue is the subtraction of the Lorcntzian 
profile for the 1+1 decays. While it is true that inte- 
grating the Lorentzian profile does give the one-photon 
formula 5Zn'~=2 -^-nl n'p f° r t ne decay rate, the Lorentzian 
profile plays no role at all in the conventional MLA. The 
derivation of the Sobolev approximation assumes that the 
line in question (here Lya) is a <5-function, and hence it 
is valid for any line profile so long as it is sufficiently nar- 
row. Thus it is not clear mathematically why in Eq. (|32p 
we should use the Lorentzian as the "already included" 
profile instead of any other function that integrates to 
Y^n~=2 Ani,n'p- A final distinct issue is that photons emit- 
ted blueward of Lya (i.e. in super-Lya decays) will be 
re-absorbed, but not instantaenously, so they do have 
an opportunity to affect recombination; we will see later 
that they yield corrections of up to ~ 1.5%. 



C. Physical resolution of the double-counting and 
re-absorption problems 

To resolve the double-counting and re-absorption prob- 
lems, in the next section we will proceed to write down 
the two-photon radiative transfer equations. A radia- 
tive transfer calculation is necessary to fully address re- 
absorption. Both numerical solution of the two-photon 
equations fSec.HV|) and analytic approximations (Sec.lV| 
will show that the two-photon decay rate and radiation 
spectrum both remain finite; the photon phase space den- 
sity f v in the optically thick region near the Lya line ap- 
proaches a quasi-equilibrium value in which the nl — > Is 
decays and Is — > nl absorptions roughly cancel. 

The double-counting problem is conceptually trickier. 
There is no physical distinction between a pure two- 
photon decay that happens to be on resonance and a 
1+1 decay, so we are free to draw the distinction how- 
ever we want so long as the total rate and the radiative 
transfer calculation are correct. Since the MLA is al- 
ready treating the 1+1 decays and associated radiative 
transfer via the Sobolev approximation for the Lyman 
lines, it makes sense for us to set up the distinction in 
such a way as to preserve the validity of the Sobolev ap- 
proximation. The assumptions underlying the Sobolev 
approximation can be identified from e.g. the derivation 
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in Rcfs. [H, |43|; they are (i) that the line width is nar- 
row enough that phase space factors, atomic density, and 
Hubble rate do not vary significantly during the passage 
of a photon through the line; (ii) the emission and ab- 
sorption profiles are the same; and (iii) that the line is 
the only significant emission or absorption mechanism in 
its frequency range. [Additionally the Lya diffusion is 
neglected, but for very optically thick lines the photons 
in the line tend to come to quasi-cquilibrium in which the 
phase space density is constant, / sa x u /(g u xi/gi — x u ). 
In this case diffusion in frequency plays little role and the 
Sobolev escape probability works well. This assumption 
will be revisited in future work where Lya diffusion is 
included.] 

Of these, (i) requires that the "1+1" decay profile to 
have support only within a frequency width Av <C Vhya ■ 
Assumption (ii) is harder: the ratio of emission to absorp- 
tion profile in a thermal medium is oc e.~ hv l T if hv ^s> T. 
Of course the real Universe is not in thermal equilibrium, 
however the emission/absorption ratio exhibits this expo- 
nential dependence in the vicinity of Lya so long as the 
excited levels n > 2 are in Boltzmann equilibrium with 
each other, even if the Saha equation is not obeyed and 
the ratio xi v jx\ s is not equal to the Boltzmann value. 
[We will prove this explicitly in Sec. IV Bl when we derive 
Eq. ([55)) : the ratio of emission/absorption profiles is f®.] 
In order for e ~ hy / T to vary by only a small fraction and 
satisfy assumption (ii), we need Av <C T/h. This "^C" 
is especially stringent: the Sobolev escape probability 
is guaranteed to be valid to accuracy e only so long as 
the emission/absorption ratio (oc e~ hu l T ) is constant to 
within e, so we need Av < eT/h. This condition may be 
relaxed if the two-photon absorption optical depth (dom- 
inated by Is — > 3d) rcdward of v\ (the red boundary of 
the Lya line) is 3> 1, since in this case all Lya photons 
that escape the line are re-absorbed anyway and the es- 
cape probability is moot. This issue will be explored 
numerically when we vary v\ . 

Assumption (iii) requires that 1+1 and pure two- 
photon decays not overlap in frequency - that is, for 
each Lyman line, we must define a frequency range 
v\ < v < V2- If a two-photon decay produces a pho- 
ton within this range, it gets counted as 1+1. Otherwise 
(i.e. if the photon is far enough from resonance) it gets 
counted as a pure two-photon decay. In this picture the 
"1+1" part of the profile resembles a Lorentzian trun- 
cated at the frequencies v± and ^2, and the far damp- 
ing wings are considered pure two-photon decays, even if 
they resemble a Lorentzian shape. 

As a practical matter there is another restriction on i>\ 
and V2- In our numerical code we are not including the 
pole displacement in the calculation of the two-photon 
transition rates, since this would introduce highly compli- 
cated temperature dependence into dk n i/dv. Therefore 
we must choose the detunings ^Lya — v\ and v-i — v-^ ya to 
be large compared to the intrinsic width r 2 p of the Lya 
line. This has the ancillary advantage that the 1+1 de- 
cay rate is not significantly reduced by the truncation of 



the damping wings: the fraction of decays that proceed 
through the far damping tails is 



/far = 1 — 



A 2p 

4^ 



l?2pdp 



^[(v - ^Lyo) 2 + (r 2 p/47r) 2 



Vhya —V\ V2 — ^ Lya 



(34) 



in the Lorentzian approximation. This quantity is <C 1 
so long as v Lya — v\ and v 2 - ^L ya are » L 2p . 

In principle the Lya line profile is affected at some 
level by two-photon corrections even within the range 
V\ < Vhya < because the decay process that emits the 
line (e.g. 3d — > 2p — > Is) interferes with other pathways 
(e.g. 3d — > 3p — > Is) that disrupt the approximation that 
there is a single pole in the matrix clement (which leads 
to a Lorentzian). This is not a problem for us: as long as 
assumptions (i)-(iii) above hold, we may use the Sobolev 
approximation for the central part of the Lya line, which 
does not depend on the line profile. 

Note that our separation into 1+1 and pure two- 
photon decays differs from CS08. Also because v\ and 
v 2 are arbitrary, one cannot uniquely define a total two- 
photon decay coefficient A„; for the two-photon decays. 
This is not a problem since the "total" A„; will not ap- 
pear in our equations and the notion of this object plays 
no role in our solution. 

Our prescription, like that of CS08, is not unique: v\ 
and V2 are arbitrary. The solution we derive can only 
be correct (and useful) if, within the restrictions above, 
we get the same answer regardless of our choice of these 
numbers. This will be verified numerically, and we will 
also find that the analytic approximations we derive have 
well-defined limits for small i^Lya — v\ and v 2 — v\^ ya that 
agree with the numerical results. 



D. Two-photon radiative transfer equations 

In the last two sections we highlighted the reasons why 
it is difficult to graft two-photon transitions onto the 
MLA. We will now write down the full set of two-photon 
equations, including those for the radiation field. This 
full set of equations must include both two-photon emis- 
sion and absorption in order to be consistent with the 
laws of thermodynamics, and should serve as a starting 
point for any numerical solution or analytic approxima- 
tion. After writing down the equations we will consider 
two ways of implementing them in an MLA. 

In the case of a two-photon decay from nl — ► Is, the 
net number of decays per unit photon frequency per unit 
time per unit H atom is 



A n l{v) 



dA n i 
dv 



(1 + /fX 1 + U')Xnl - —fufu'Xu 

9u 



(35) 
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so that the net two-photon decay rate is 

n(l-n~ 2 )TZ 

iniU-, = - / A„/(i/) dv 

J(l-n- 2 )TC/2 



for n > 2, and 



n£,n>2 



(36) 



(37) 



In the case of the 2s — > Is transition, these equations 
replace Eq. (fTTj) . Note that it is only necessary to fol- 
low the higher-energy photon (y > is'): the lower-energy 
photon has v' < 1Z/2, and hence has negligible proba- 
bility of re-exciting a hydrogen atom. We can calculate 
the probability of these photons interacting again by con- 
sidering their optical depths for various processes. The 
optical depth per Hubble time for two-photon absorp- 
tion is typically in the range 10 -20 -10 -16 at frequency 
= 3TZ/8. The optical depth per Hubble time for 
Balmer continuum absorption is higher, reaching a peak 
of 2 x 10~ 4 at the threshold 11/4 and z = 1400. The 
optical depth in the Ha line is 5.5 x 10 -4 at z = 1400, 
and since this is a resonance in the Raman scattering 
cross section the optical depth for Raman scattering is 
much less (typically by a factor of the resonance width 
divided by the detuning from resonance). Since there 
are ~ 0.5 net two-photon decays per hydrogen atom dur- 
ing H recombination, the latter numbers suggest that 
re-absorption of the lower-energy photon from the two- 
photon decays would influence recombination at the level 
of a few parts in 10 4 . In practice the effect should be less 
since at z ~ 1400 the excited states H(n > 2) and the 
continuum H + + e~ arc kept in equilibrium (in which case 
absorption or emission of Balmer photons has no effect) 
and at lower redshift the lower population of the 2s and 
2p levels results in much smaller optical depths. There- 
fore we have not followed the lower-energy (i/) photons. 

Because some portions of the two-photon spectrum are 
optically thick, we will also need the equations for the 
photons. For frequencies not on an H I resonance line, 
the photons follow the equation 



U=Hv 



c 3 n H 



df 8TTV 2 



]Ta„^), 



(38) 



where the first term is the Hubble redshift, and the sec- 
ond term is the source (the coefficient describes the pho- 
ton phase-space density). 

One can write similar equations for two-photon recom- 
bination/ionization and Raman scattering. Only those 
reactions involving the ground state are included because 
one can transition among the other states by optically 
thin one-photon transitions, which are much faster. In 
the case of the Raman scattering reaction: 



we extend Eq. (|35p by writing 
dK nl 



A n l{v) 



dv 



(1 + fv)fv>Xnl 



9nl 
9U 



+ /„')*! 



(40) 

for v > (1 - tT 2 )TI and v' = v - (1 - tT 2 )TI. Here 
dK„i/dv is the differential Raman scattering rate, de- 
fined by Eq. (|B10[) . Equation (f3"6"| is extended to include 
Raman scattering by writing 



•Enl | Raman 



l-n- 2 )K 



A n i(u) dv; 



(41) 



we have only extended the integration up to 1Z because 
if the output photon has a frequency above 1Z then it 
immediately ionizes an H atom. For two-photon recom- 
bination to the ground state, i.e. 



H+ + e" -> H(ls) +hv + hv', 



(42) 



one writes 



A c (v) 



(1 + /„)(! + f„>)nH_x e x p PM(E e ) 



—f v f v >— ]aW(v,v')dE e , (43) 

dEdV 2 J v ' ' ' v ' 

where hv 1 = E e — E\ s — hv, and a^(v, v') is the differ- 
ential two-photon recombination rate coefficient (in e.g. 
cm 3 s _1 Hz -1 ), given by Eq. (|B15|) . The A c (v) term is 
added along with the bound-bound two-photon transi- 
tions A n i (v) in Eq. 



IV. NUMERICAL APPROACH 



H(nZ) + hv' -> H(ls) 



hv, 



(39) 



In Section IIII Dl we wrote down the radiative trans- 
fer and level population equations including two-photon 
transitions. Now we solve them numerically by discretiz- 
ing the two-photon continuum. We first write down the 
methods of solution fSection llV Ap and describe our fidu- 
cial choice of parameters (Section IIVBI) . Then we con- 
sider the results obtained by this method in Sect ion HV CI 
turning on a sequence of effects in order: (A) the "cor- 
rect" handling of the 2s — -> Is decay, including stimu- 
lated emission and feedback; (B) sub-Lya two-photon 
decays, i.e. those from n > 3 where both photon fre- 
quencies are below Lya; (C) super-Lya two-photon de- 
cays, i.e. those where one photon frequency is above 
Lya; (D) Raman scattering; and (E) two-photon recom- 
bination/ionization. 



A. Method 

Our agenda here is to discretize the two-photon trans- 
fer equations. We will do this by converting the continu- 
ous frequency distributions such as dA n i/dv into discrete 
distributions, i.e. lines. We will define a set of virtual 
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levels, such that a two-photon decay such as 3d — > Is 
can be represented formally as a one-photon decay from 
3d to a virtual level, and then from the virtual level to Is. 
We will then write down a set of equations for these lev- 
els so that the equations for the occupation of the virtual 
level are formally equivalent to the two-photon radiative 
transfer equations. We emphasize that this is only a for- 
mal correspondence used to force an MLA code to solve 
the two-photon transfer equations, and that the virtual 
levels are not actual levels of the hydrogen atom (they 
are not linearly independent states in the Hilbert space) . 
In order to be pedagogical, we will work through only 
the two-photon decay/excitation here, and then state the 
analogous results for Raman scattering and two-photon 
recombination / ionization. 

In order for the discretization procedure to work, the 
total rate and (after appropriate smoothing) shape of 
the two-photon distribution must remain fixed. Since 
two-photon absorption also plays a key role, we must 
also make sure that the two-photon optical depths 
and branching probabilities (relative probability of two- 
photon absorption to different excited levels) are pre- 
served. 

The simplest way of doing this is to multiply each prob- 
ability distribution by a discrctizing function D'\ 



dA r , 



dv 



dk„i 



used 



dv 



D'(u), 



(44) 



with 



N v 
b=l 



Vb), 



(45) 



with the "spike" frequencies ordered v\ < v^ < ... < 
z//v„ and with Av b being the spacing between consecutive 
spikes. One can also define the integral 



D{u) 



D'(v) dv w v — v c 



(46) 



As always the delta function has infinitesimal width. 
The frequency v cut is the minimum frequency at which 
we follow the spectral distortion. Note that dA n i/dv is 
only defined above half the maximum frequency, v > 
(1 — n~ 2 )7Z/2; below this cutoff we will formally take 
dAni/dv = 0. We take /„< to be a blackbody here, i.e. 



we neglect the spectral distortion at v' < (1 



l )TZ/2. 



As shown in Section fill D[ this does not lead to any sig- 
nificant errors. 

Now that we have a discretized spectrum, we proceed 
to solve the radiative transfer equations. At frequencies 
in between spikes, Eq. (|38p is simply a free streaming 
equation and has a solution analogous to Eq. (|13[) : 



f(v b + e,z) = f 



Vb + 1 / 1 i \ i 
Vb+i - e, (1 + z) - 1 

Vb 



(47) 



Within a spike the situation is more subtle. At v w v bl 
Eq. ([35]) becomes 



fu 



Hv b — + o 2 D \ V ) 

dv 8irv£ 

nl 



dv 



E 

7 U 



9nl 
Sis 



fv'Xu - (1 + U')x n l 



dA 



III 



dv 



Since the low-frequency photon comes from the CMB (i.e. 
we are neglecting two-photon absorption where the low- 
energy photon is part of the spectral distortion), we may 
take /„< to be constant across the spike (but of course it 
does depend on n). Then we may define the coefficients 



9nl 
Sis 



fv'Xu - (1 + f,y)x n i 



dA nl 
dv 



and 



Bo =$> + /"') 



dA r , 



■n I 



dv 



Then Eq. 

K 

Hv b 



can be written as 



df v 



c 3 n H 



dv 8irHv3 



D'(y){B -B 1 f v 



(49) 



(50) 



(51) 



This equation has an integral solution, 
'Bo 



Si 



-exp^^ s B 1 [D(y)-D(u b 



h 

Hv b 



\-£^B x [D{y)-D{y)] 



(52) 



As we make the delta functions infinitesimally narrow, 
the integrand in the last term remains finite but the range 
of integrated variable becomes infinitesimal. Therefore 
this term disappears. We can simplify Eq. ([52")) by defin- 
ing 



Tb{v) = - 



8irHv3 



Bi[D(i/) - D(yb + e)}, (53) 



and 



Ar fc = T b (v b - e) = 



8-kHvI 



BiAvb. 



(54) 



Comparing with Eq. |49|) this looks similar to a Sobolev 
optical depth; indeed, Av b is the optical depth to two- 
photon absorption. 
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At the end of the day the relevant numbers to know 
are the net two-photon decay rate from each excited level 
via bin b, x n i\2-y,b, and the photon phase space density at 
Vb — (■ The latter can be computed from Eq. (|52p ; it is 
given by 



fu 



(55) 



The two-photon decay rate can be written as (see Eq. [35 



(1 + /*/)(! + fu')Xrd - —fvfv'XU 

9u 



dA„i . 
x— — D (y)dv. 
dv 



(56) 



Using Eq. (|52|) for the phase space density within the 
spike, and recalling that the last term is zero, we see 
that 



Xnl\2~(,b 



dv 



[(1 + + fu')x n i - f Vb f v >Xi a )Av b , 

(57) 

where the weighted mean phase space density has been 
defined as 



f -* + (f B °\ 1 ~ e ~ AT " 



(58) 



It is convenient to define the Sobolev-like "escape prob- 
ability" from the spike, 



i 



-Ar 6 



A n 



(59) 



In principle we could directly add Eqs. (|55p and 
Eq. (|57|) to the level equations, since f Vb is linear in 
the excited state populations so long as they are small 
enough that x\ s dominates the sum in B\ (Eq. [49]) . We 
want to do this in a way that is easily incorporated into 
the steady-state equations; this can be done by invent- 
ing a "virtual" level b for each spike, and extending the 
steady-state source vector s and transition matrix T to 
include it in such a way as to reproduce the two-photon 
equations. The populations of the virtual level for each 
spike will correspond directly to the spectral distortion, 
which will be convenient numerically. Let us define for 
each spike an arbitrary "degeneracy" gb- Formally one 
should think of gt as infinitesimal, since the virtual levels 
cannot be occupied, but within the context of the steady- 
state approximation the argument below will work for 
any arbitrary value. We will then define the virtual level 
occupation 



_ 9± 7 
Xb — 2 x lsju b 



and the virtual level energy 
Eb = Ei s 



hv. 



(60) 



(61) 



We will also define the virtual rate coefficient 

dA„,i 



A, 



Lb 



dv 



Av b . 



(62) 



Then Eq. (|57|) gives us the effect of a given spike on the 
rate equation for level nl. Since /„< = f n l,b+ an( i fu b <C lj 
this is 

X n l\b = -A n i b {\ + f n l,b+)x n l + —Anibfnl,b+Xb- (63) 

9b 

This is exactly equal to what Eq. (|18|) for Rji would give 
when plugged into the steady state equation, Eq. p^]) . 

In order to implement something like the steady state 
equations, we must also express the equation for f Vb 
(Eq. I5"8")) in language that allows us to define a T-matrix. 
That is, if we imagine extending Eq. (|2~Tj) to include both 
real (r) and virtual (v) levels, 



rj-irr rprt) 



(64) 



we want to construct values of s v , T vr , and T ™ that are 
equivalent to Eq. (|5"5|) . (Note that we have already con- 
structed the top row of the matrices.) Re-writing Eq. (|58|) 
in the form 



9nl 

9b 



Anlh X b 



\ nl 

+(i-n 6 )^(i + /^)^ 

i il 

3b J gu 
we can then define the virtual level source term 



fs^f 9nl a i 9b f n fc 
VtT 9b ) gis i - n b 



the virtual-real sub-block of the T-matrix, 

Tb.nl = — + fu')Anl,b, 

nl 

and the virtual-virtual sub-block, 

<W £ 9m 



Tbb> — 



l-U b 



^ 9b 



Lb- 



(65) 



(66) 



(67) 



(68) 



nl 



A simple algebraic calculation shows that this is precisely 
what one would calculate for sj,, Tb >n i, and Tw if one were 
to make the replacement 



AbAsPbSs 



n fc \ - g n i , . 

Z ^j- / Jv'A n l b 

i - n 6 ^ g b 



(69) 



in Eqs. (fl~7|) and ((20|) . This can also be re- written as 

1 - e- Arb 



AbAsPbAi 



c 3 n H 1 - (1 - e- AT »)/An' 



(70) 
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In other words, the steady-state equations can be ex- 
tended to include the two-photon decays by adding a set 
of virtual levels, inserting a virtual one-photon transition 
between each level nl and a virtual level with Einstein co- 
efficient given by the two-photon decay rate, and intro- 
ducing a transition from the virtual level to the ground 
state Is with rate given by Eq. ([TO)) . 

The above methodology can be extended to in- 
clude Raman scattering and two-photon rccombina- 
tion/ionization with no conceptual changes but with the 
addition of some new reaction rates. The Raman scat- 
tering reaction, Eq. (|B7[) , is included by including a new 
radiative rate 



■ 9ni dK nl 

Ab.nl = ; — Az/fc, 



!-)b 



dv 



(71) 



and two-photon recombination is included by writing an 
effective recombination coefficient to the virtual level: 



a b (E e ) = a {2 \v, v')Av b , 



(72) 



where hv' = E e — E\ s — hv. Equation (|70[) remains valid 
as long as we define a total optical depth 



a c 3 n H [v-^ 9nl , 



III 



dv 



dK n , 
dv 



+ 



o7/2 3/2 poo p -, 

' dE e aW(v,v')^ f u ,dE e (73) 



9u 



B. Parameters 

We have run our numerical radiative transfer code 
several times with different sets of parameters. These 
choices are summarized in Table [I] The truncation fre- 
quencies of the Lya line (outside of which two-photon 
decays are treated using the numerical radiative transfer 
instead of being treated as 1+1 decays with the Sobolcv 
approximation) are v\2 = v-^ ya ± Az/ max . In the lat- 
ter case one should compare to the intrinsic width of 
the 2p intermediate state, I^p = A2 P ^ S = 6.3 x 10 8 s _1 
in vacuum. (The additional finite-temperature broad- 
ening is due mainly to 2p — ► 3s, 3d transitions and for 
recombinaton-era temperatures is -C ^-2p,is-) Thermal 
equilibrium conditions were assumed for the initial data 
at z; n it. This redshift was chosen to be after the com- 
pletion of He II— >He 1 recombination in order to avoid 
needing to simultaneously follow hydrogen and helium. 
The "Basic" set of parameters will be used for all plots 
in the paper unless otherwise indicated. All runs com- 
menced at Zmit = 1605.8 and used n max = 30 shells of 
the hydrogen atom (465 sublevels including I resolution) . 

For the Basic run, the frequency grid extends only up 
to Lyi = 0.997?. since above this frequency the real hydro- 
gen lines are very closely spaced and feedback is almost 



TABLE I: The parameters for each run of the radiative trans- 
fer code. Shown are the number of virtual levels iVvirt; and 
the step size in log scale factor A In a; the maximum detuning 
from Lya at which a photon is considered resonant, Av max ; 
the step size in the bound-free rate integrals, Aln_E e ; and 
the line corresponding to the maximum frequency used in the 
grid. 



Run 


iV v i rt 


10 5 Aln 


a Ai/ max 


A In E e 


Max. grid 








GHz 




frequency 


Basic 


367 


4.25 


105 


0.1 


Lyt 


WideLinel 


365 


4.25 


315 


0.1 


Lyt 


WideLine2 


365 


4.25 


945 


0.1 


Lyt 


LoRes 


262 


4.25 


6116 a 


0.1 


Lyt 


HalfStep 


367 


2.125 


105 


0.1 


Lyt 


HiRes 


695 


4.25 


105 


0.1 


LyM 


CoarseBF 


367 


4.25 


105 


0.2 


Lyt 


"For the red 


boundary v\ ; 


the blue 


boundary 


is at I'Lya 



10259 GHz 



instantaneous. The "WideLinel/2" runs are the same 
as the Basic run except that we increased Ai/ max by a 
factor of 3 and 9 respectively, i.e. a two-photon decay 
had to emit a photon within 315 (instead of 105) GHz 
of the Lya line in order to be considered a 1+1 decay 
instead of being treated with the two-photon radiative 
transfer. Its purpose is to show that Az> max does not 
matter. The "HalfStep" run is a convergence test that 
used half the time step of the Basic run. The "LoRes" 
run has a lower-resolution grid of virtual levels, while 
"HiRes" has a higher-resolution grid (twice the resolu- 
tion except within 25 THz of Lya) that extends all the 
way up to Lyn (Is — 13p). The "CoarseBF" run uses a 
step size of A\nE e = 0.2 instead of 0.1 for evaluating 
the energy integrals for the bound-free rates. All of the 
alternative runs gave changes in the ionization fraction 
\Ax e \/x e < 10~ 4 , except WideLine2 and LoRes for which 
the maximum change was 4 x 10~ 4 . 



C. Results 

The change in the recombination history caused by 
the new physics is shown in Figure [3J We will discuss 
the effects qualitatively here; analytic estimates will be 
given in Section [V] and compared with the numerical re- 
sults. The strength of the numerical result is that the 
only approximations made are those intrinsic to the two- 
photon radiative transfer equation and the step sizes, fre- 
quency resolution, etc., which can be varied to establish 
convergence. The analytic approximations yield insight 
but require numerous additional simplifying assumptions 
whose validity is difficult to assess. 

The first effect that we have considered is the inclu- 
sion of stimulated two-photon decay and nonthermal two- 
photon absorption in the 2s — > Is transition (A). These 
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FIG. 3: Upper panel: the recombination history. Middle 
panel: The absolute change Ax e caused by incorporation of 
each process. Lower panel: The relative change Ax e /x e . The 
changes are cumulative in the sense that the "B" curve rep- 
resents the effect of incorporating A and B, etc. The (D) and 
(E) curves are almost indistinguishable on the plot. 



effects go in opposite directions since stimulated emis- 
sion increases the effective 2s — ► Is decay rate and hence 
accelerates recombination [23 whereas nonthermal two- 
photon absorption retards it [24[ . Initially the net effect 
is a slight acceleration because of the strong stimulated 
emission at high temperature, but by z = 1345 this turns 
into a retardation (Ax e > 0) as enough hydrogen recom- 
bines to build a significant Lya spectral distortion. This 
grows to Ax e w 0.001 at z w 1200 and then tapers off. 
For most redshifts we are in agreement with Ref. pij 
(compare their Figure 9 to our Figure [3]), but at the tail 
end of recombination (z < 800) we find a smaller effect, 
Ax e /x e = 0.0018 at z = 700 instead of 0.0047. This is 
probably due to the approximation in Ref. [24| that the 
excited levels n > 2 are in Saha equilibrium; if we force 
Saha equilibrium in our code, our Ax e /x e rises to 0.0048. 

Next we considered (B), the sub-Lya two-photon tran- 
sitions ns,nd <-> Is. In these cases the emitted photons 
can never be resonantly absorbed by a hydrogen atom. 
As one would intuitively expect, the addition of this new 
process brings the Universe closer to thermal equilibrium, 
that is it accelerates recombination. These transitions 
have a huge effect: they reduce the electron abundance 
by as much as ~ 2%, and moreover this reduction is oc- 
curring at z ~ 1100, i.e. at the surface of last scattering 
where it matters most. 

The next effect is (C), the super-Lya two-photon de- 
cays. This is another process connecting excited states 
of H I to the ground state and naively we would expect 
that like (B) it would accelerate recombination. This 
is true at first. However at z = 1290 the situation re- 
verses and the new process delays recombination. The 
reason for this is that these two-photon decays are pop- 
ulating the region of photon phase space at v > VLya 
to far above the blackbody value. In the later stages of 
recombination, these photons redshift down to Lya and 
begin re-exciting hydrogen atoms. This delayed feedback 
is present even in the standard recombination calculation 
because some photons escape from the Ly/3 line [4(| but 
the two-photon decays to frequencies between Lya and 
Ly/3 make a significant addition because their probabil- 
ity of immediate re-absorption is much lower. A similar 
result occurs for (D) Raman scattering, which initially 
accelerates recombination but then slows it down. In this 
case the main effect is 2s — > Is scattering in which the 
incoming photon is at v < and the outgoing photon 
is at is Lya <v< VLyp. 

Two-photon recombinations (E) have an accelerating 
effect as one would intuitively expect. However since the 
electron occupation probability of the continuum states 
(E e > 0) is much less than the 2s, 3s, and 3c? states this 
is a small effect; we find a maximum of |Aa; e |/x e ~ 10~ 3 . 

We can also examine the effect on the high-frequency 
end of the spectral distortion. In Fig.QJ we show the pho- 
ton spectrum at three epochs during recombination. The 
full calculation inlcuding effects A-E is shown compared 
to the case with effect A (modifications to 2s *-* Is) only. 
(We show the case with the 2s *-* Is modifications turned 
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FIG. 4: The radiation spectrum as a function of time. Each panel shows the full spectrum (solid line); the spectrum neglecting 
the highly excited two-photon and Raman transitions (i.e. including only modification A) (dashed line); and the blackbody 
spectrum (dotted line). We do not show the lower frequencies because our code does not follow the lower-energy photon in 
two-photon decays nor the Balmer, Paschen, etc. lines or bound-free continuum. 



on because in the original case with the 2s — > Is decays 
followed with the integrated rate coefficient A2 S only, the 
code cannot follow the two-photon spectral distortion.) 
The large jumps at each Lyman line are due to photons 
redshifting out of the line, as described by Peebles @. 
However even at early times there is a spectral distortion 
between the Lyman lines that is due to two-photon and 
Raman transitions. Particularly impressive is the phase 
space between Lyev and Ly/3: here the photon occupation 
number can be up to an order of magnitude larger than in 
the traditional treatment. Due to their very long mean 
free paths, these photons can survive to z < 1100 and 
slow down recombination by re-exciting hydrogen atoms. 



V. ANALYTIC APPROACH 

In the previous section we numerically evaluated the 
effect of two-photon transitions on hydrogen recombina- 
tion using a radiative transfer scheme based on "virtual" 
levels. Here we consider simple analytic estimates of their 
effect. The main purposes are to provide an indepen- 
dent check of the fully numerical result and to develop 
some understanding of the magnitude of the various ef- 
fects. We consider the improvements to the 2s <-> Is 
decays (Section IV A[) ; sub-Lya two- photon decays (Sec- 
tion |VB|; super-Lya two- photon decays (Section |VC[) ; 
and Raman scattering (Section lVDp . We also considered 
the two-photon recombinations in the numerical code, 
however since they were found to be unimportant we will 
not construct an analytic estimate. 

Our analytic investigation of the two-photon processes 
is based on many simplifying assumptions, most notably 
that (1) T r -c MZ; (2) the low-lying excited states of 
the hydrogen atom (n = 2, 3) are in Boltzmann equilib- 
rium; and (3) that the radiation field can be treated as 
being in steady state (or some lowest-order non-steady 
state approximation, in the cases of two-photon decays 



where one photon is above Lycv and Raman scattering). 
These assumptions are only marginally satisfied during 
recombination and it is difficult to directly estimate the 
residual error. None of these approximations is made 
in the numerical code. For these reason we believe the 
numerical results in the previous section are the most 
reliable estimate of the effect of two-photon transitions. 



A. Stimulated emission and nonthermal absorption 

in 2s <-» Is 

The 2s *-y Is transition is modified by both stimu- 
lated emission and absorption of the Lyman-a spectral 
distortion. The approach described here, like that of 
Refs. [23|, [24[, assumes the optically thin limit for the 
Is — > 2s absorption, but also takes the leading-order 
term in the power series for the spontaneous decay rate 
in order to arrive at an analytic formula. This effect re- 
quires a thermal CMB photon to be present at one of 
the photon frequencies v' and hence occurs in the regime 
hv' ~ T r -C MZ. In this regime the 2p intermediate state 
dominates; using the exact hydrogenic matrix elements 
(ls||r||2p) = 2 15 / 2 a /3 9 / 2 and (2a||r||2p) = -3 3 / 2 a , we 
find (sec Eq. [B4| 

dMs „ 512a| v> 
dv ~ 729 K' K ' 

The total rate A:fc| of stimulated emission minus ab- 
sorption of the spectral distortion is 

... /-fLyo/2 J A 

Ai l = J q f,' (x 2s -x ls A f u )dp', (75) 

where A/„ is the phase space density of spectral distor- 
tion photons (i.e. the total phase space density with the 
blackbody subtracted) . Since /„< is a blackbody function, 
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we may write 



A.i 



. (A) _ 512og, f^ /2 v'(x 2s - x u Af v ) 



729ft 



e hv'/T T 



dv'. (76) 



Here x 2s and x\ s are independent of frequency, and Af v 
varies smoothly with v' (since we are looking on the red 
side of Lya; if we considered negative v' there would be a 
jump) . Therefore this integral can be approximated using 
one-point Gaussian integration (see e.g. Chapter 19 of 
Ref. [il] .) The idea of one-point Gaussian integration is 
that if one has an analytically integrable weight function 
w{v), and a linear function g(y') = a + bv', that: 



w(v')g(v')dx = C g(v' ), 



(77) 



where Cn = / w(v')dv' and v = Cq 1 j v'w(v')dv' . This 
enables integration with a single sample of the function 
g instead of two as is usual for linear functions. For 
functions g{v') that vary smoothly and are being inte- 
grated over a narrow range we may use Eq. (|77|) as an 
approximation. In our case, w{v') = v' /{e hu l Tr — 1) and 
g{v') = X2s — x\ s Af v . Application of Gaussian integra- 
tion gives: 



L 



V ^«' 2 v'{x2s-X ls &f v )dv' 



in 
where 



and 



ohv' /T r _ i 



C (x 2s - x ls Af„ ), 



Co 



v'dv' 



n 2 T 2 



Q e hu'/T t _ 1 Q h 2 



00 - ,2 dv' 



Q e hv>/T r _ l 



12C(3)T r 
ir 2 h ' 



(78) 
(79) 

(80) 



(The upper cutoff in the integral does not matter because 
in practice the exponential provides a cutoff.) Here £(3) 
is the Riemann (^-function, and vq = vi^ya — v'q. This 
leads to 



Ai (A) ^ 256^ 2 a|T r 2 



2187/i 2 ^ 



(x 2s -x ls Af Uo ). (81) 



We have compared Eq. (fST|) to the full numerical result 
for the 2s <-> Is two-photon transitions in Fig. [5] 



B. Sub-Lya two-photon decays 

We now consider the sub-Lya two-photon decays from 
higher excited states of H I such as 3s and 3d. In most 
cases the higher-frequency photon is just below Lya and 
can be cosidered to be part of the red damping tail. 
Therefore in order to construct an analytic theory of the 
two-photon transitions we must develop an approxima- 
tion for the radiative transfer in the red damping tail. 



A. Effect from modifications of two-photon decays from 2s 
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FIG. 5: A comparison of methods for handling the 2s «-> 
Is two-photon transitions. We show the numerical radiative 
transfer method (solid line) and the analytic approximation 
of Eq. (J8TJ) (dashed line). 



The approximation described here draws heavily on that 
of Appendix B of Hirata & Switzer ■ 

The two-photon transfer equation, Eq. with 
Eq. (|35p used for A n i (y) , can be written as a linear equa- 
tion for /: 



Hv ~ dv KvKU lvh 
where we have defined the functions 



c u ?i H 



8ttHis 3 t-f 

nl 



— fu'Xis - (1 + f v ')x n l 

gu 



dA nt 
dv 



and 



f° 



Snt( 1 + /^) a: "' 

T,nl[(9nl/9n)fu'Xl s ~ (1 + fu')x n l] 



(82) 



(83) 



(84) 



Here k v can be thought of as an optical depth per unit fre- 
quency (this is more convenient in cosmology than optical 
depth per gcm~ 2 as is usual in stellar structure), and /° 
is the photon phase space density that one would have 
if the level populations were maintained but the density 
were taken to oo (high optical depth limit). Note that 
the "optical depth" used here is only the optical depth 
to two-photon transitions and does not include Rayleigh 
or Thomson scattering. 

Our goal here is to develop an analytic approximation 
to Eq. (|82|) valid in the red damping wing of Lya. To do 
this we make three assumptions. We begin by making 
the steady-state approximation that the time derivative 
on the left-hand side can be dropped. Physically this cor- 
responds to the statement that the radiation spectrum in 
the vicinity of Lya changes slowly compared to the time a 
photon takes to redshift out of the optically thick part of 
the damping wing. Secondly, we assume that the excited 
levels in > 2) are in equilibrium with 2p, i.e. x n ijxi v is 
given by the Boltzmann factor [(2l + l)/3]e^^ Enl ^ E2 "^ Tr . 
This, combined with x 2p <SC x\ s , allows us to write 



ft 



X2p h(v- 



^Lyc,)/T r 



3xi 



(85) 
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Our third assumption is that we may approximate 
dA n i/dv by the single pole approximation, dk n i/dv oc 
(y — t'Lya) -2 , in which the two-photon differential decay 
rate is proportional to the inverse-square of the detun- 
ing from the Lya resonance. Because of the pole in the 
two-photon matrix element M. from the 2p intermediate 
state, this is the leading-order term in the power series 
expansion of dA n i/dv around v = v^ ya . The coefficient 
of this pole is taken from Eq. (|B5[) : in the vicinity of Lya, 
we have 



dA n 



27a?u 



dv CLqR. 6 {v — ^Lya) 2 

It follows that 

w 



|M||r||2p)(2p||r||l S )| 2 . (86) 



- v Lya y 



(87) 



where 
W 



E 



21 



c 3 nnxu 

Z P^ S nl,n>3 



|H|r||2p)(2p||r||ls)| 2 .(88) 



[We have neglected x n i in Eq. (|83|) since even with the fac- 
tor of 1 + /„/ instead of f v i it will be negligible compared 
to X\ s ; this is physically equivalent to neglecting two- 
photon decays where the higher-energy photon is stimu- 
lated.] This expression can be reconstructed in terms of 
Einstein coefficients: 



W 



c 3 nnA 2 p^ s xi s 



E ■ 

n/,n>3 



21 + 1 



A 



nl,2p 



(89) 



The prefactor is recognizable as the Lya (Gunn- 
Peterson) optical depth, 



T Lya 



_ 3c 3 n B _A 2 p^ s xi s 



so we write 
W 



nl,n>3 



21 + 1 A nl . 



2p 



ohv n t,a p /T x _ i ' 



(90) 



(91) 



With this set of approximations, the radiative transfer 
equation reduces to 



(1U = w 

dv (P - Vhya) 2 



M"-^yc.)/T t 



X2p 



3.c 



Is 



(92) 



The change of variables y = h(v — fL ya )/T T and /„ 
X2 P $(y)/(3xi s ) gives 



dy 



W 



(93) 
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FIG. 6: Solid line: The correction to the Lya escape proba- 
bility, $(— oo), as a function of W . Dashed line: The integral 
X that describes the number of photons on the blue side of 
Lya emitted in two-photon decays. Note that as W — » both 
$(-oo) — 1 and I — > 0. 



where W = hW/T r is a dimensionless number that quan- 
tifies the normalization of the two-photon opacity. We 
can understand its physical significance from Eq. (|87|) . 
For W <C 1 the integrated two-photon optical depth at 
frequencies less than some critical value v c < fLya is 



k v dv 



W 



h VL ya 



(94) 



so two-photon absorption is optically thick out to a de- 
tuning WT r /h from the Lya line and in the optically 
thick region the exponential factor in Eq. (|92|) is a pertur- 
bation. The physical situation realized in recombination 
is always W < 1 - we find a maximum of W = 0.034 - but 
the corrections due to non-negligible W are important for 
the accuracy desired in upcoming CMB experiments. 

Equation (|9"3"|) forces the boundary condition $(?/ = 
0) = 1 , and the equation can be evolved by integrating to 
the left (negative y) . It is readily seen that as W — > this 
boundary condition at <J>(y = 0) forces $(— oo) — > 1. The 
numerical integration of Eq. (|93j) with an implicit ODE 
solver is extremely fast compared to the rate coefficient 
evaluations and linear algebra required when solving re- 
combination, so we have simply solved Eq. (|93[) each time 
a value of <&(— oo) is needed. [Since $(— oo) depends only 
on the single independent variable W one could generate 
a table of values and use spline interpolation if the need 
arose for a faster evaluation.] The numerical solution is 
shown in Fig. [5] 

The net rate of two-photon transitions to Is that in- 
volve photons in the red tail of Lya can be obtained 
as follows. We compare Eq. (f3"6"|) , which gives the net 
nl — > Is rate via these transitions, to Eq. (f38|) . which 
gives the rate of production of photons in the red tail of 
Lya. Integrating the latter from frequency v\ to fh ya , 
and making the steady state approximation again to re- 
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move f v , gives 



= H 



v—dv + > 
av ^-^ 



ni 



V ^ a C 3 71 H 



A nl (v)dv. (95) 



Now if almost all of the two-photon decays from the n > 3 
levels emit one of the photons just below Lya, we may 
take v « Vhya in the prefactor of the first integral and 
the denominator of the second integral. Then 



3 

= Hv Lya (f^ ya - f Vl ) + - — ^- V / A n i(v)dv. 



E 



8™Lya 

(96) 

The last integral is the net rate of two-photon decays 
with the higher-energy photon emitted between v\ and 
Lya. (Here "net" means with the two-photon absorp- 
tions subtracted.) Since A n i(v) is strongly peaked near 
Lya due to the shape of the two-photon spectrum, we 
may take v\ to be redward of the significant two-photon 
emission, equate the last integral with £ni|27,^<i; Ly „, and 
then have 



= Hv\, ya .(jv 



l T3 \ 

where Ax^ is the net rate of two-photon transitions to 
Is that involve photons in the red tail of Lya. Recalling 
the definition of &(y), that ^L yQ and v\ correspond to 
y = and y = — oo respectively, and that $(0) = 1, we 
find 

a %2p 



c 3 n H a -(B) 



(97) 



Ax 



(B) 

1 



c 3 nn 



3xi 



[$(-oo) -1]. 



(98) 



This analytical approximation is not a surprise: it is just 
an extension of the Peebles [f| argument that the net de- 
cay rate via emission/absorption of photons in a narrow 
region of the electromagnetic spectrum that is taken to be 
in steady state is equal to the difference in phase space 
densities on the two sides of the region, multiplied by 
the appropriate normalization factors. In fact, Eq. (|98p 
without the factor $(— oo) — 1 (and with a correction 
for blackbody photons entering on the blue side) is the 
Peebles formula for the Lya decay rate, so we can think 
of $(— oo) as representing an effective enhancement of 
the Lya production rate due to two-photon transitions 
into the red damping wing. In our approximation this 
enhancement factor depends only on W. 

We now have all of the pieces required to incorporate 
the additional decay rate given by Eq. (f9"8")l into our re- 
combination model. Rewriting the prefactor of Eq. ([98]) 
in terms of TLya, we can simplify Eq. (|98[) to 

An 



Ax 



(B) 



T Lya 



£2p[$(- 



(99) 



where $(— oo) is evaluated as a function of W. Note 
that in the W — » limit, corresponding to negligible two- 
photon opacity, $(— oo) — > 1 and this correction vanishes. 
We have compared Eq. (|9"9"|) to the full numerical result 
for the n > 3 two-photon transitions with v < fj_,y a in 

Fig.m 
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FIG. 7: A comparison of methods for handling the sub-Lya 
n > 3 two-photon decays. We show the numerical radiative 
transfer method (solid line) and the analytic approximation 
of Eq. (99} (dashed line). 



C. Super-Lya two-photon decays 



Wc now turn to the super- Lya two-photon decays from 
n > 3 levels of H i. These decays arc fundamentally 
different from the sub-Lya decays considered in the pre- 
vious section. The reason is that a photon emitted at 
frequencies below Lya has a nonzero probability of es- 
caping the Lya line, after which it is extremely unlikely 
to be re-absorbed; it may be Thomson or Raylcigh scat- 
tered, but it is unlikely to re-excite an H I atom and 
will consequently have little further impact on recombi- 
nation. In contrast, all photons emitted above the Lya 
frequency will eventually redshift down to Lya and be 
absorbed. That is, for every super-Lya two-photon de- 
cay there is an additional Lya absorption, so the implied 
net number of decays to the H I Is state is zero. For 
this reason, CS08 did not include these decays in their 
analysis. However, we can see from Fig. [3] that these de- 
cays do indeed result in a change in the recombination 
history. This is because of the delay between the two- 
photon decay that produces the photon at v > v-^ ya and 
the Lya absorption. This section constructs an analytic 
approximation that accounts for the effect of this delay 
on recombination. 



, 2 7 



The key to our analytic approximation is to consider 
, the number of non-thermal photons above the Lya 
frequency per H atom produced by the two-photon de- 
cays. Initially the photon spectrum is thermal, so x 2 ^ = 
0. Long after recombination is over and all of the photons 
have redshifted below Lya, we also have x 2 ^ = 0. The 
net rate of decays to the ground state due to two-photon 
decays with v > vx^ ya and re-absorption of the emitted 
photons is 



AC) 



(100) 



19 



The total number of such decays per H atom is 

;,( C ) r H - I ~2 7 M _ J2 7 ,. _ 2-y 



\ dt = - J ± + 7 dt = X + 7 (ifi„al) - £ + 7 (tinit) = 0, 

(101) 

as expected. 

It only remains to find x2(t). By mode counting, we 
find 



C 3 ?lfJ 



(102) 



where A/„ is the non-thermal phase space density of 
photons. Switching to the variable y = h{y — VLy a )/T T , 
and taking the approximation that the distortion photons 
above Lya are mostly in the regime v/vi^ya — 1 <C 1, 

^ 8tti/ t 2 v T r f°° 



We may estimate the spectral distortion using Eq. ([82 
Rewriting it in terms of the distortion A/„, we find 



Hv di 



- K v (Af u - A/°), (104) 



where 



(105) 



and the optical depth per unit frequency n v is the same 
as that from Scc. lVBl Under the same assumptions from 
Sec. |VB] that led to Eq. (85]), we find 



A/ o w (X3p_ _ e -h^/%\ e -h(u-^ ya )/T^ (1Q6) 



We now assume the distortion is in steady state, i.e. that 
the left-hand side of Eq. (|104| is zero. Defining 



*W = ( ^ - e-^ L - /T ' ) A/„, (107) 



we convert Eq. (| 104[) into 

0= 5-- e^). 

Using Eq. (|57|) for k„ we convert this into 
(7* 



dy 



W(e v ^ - 1). 



(108) 



(109) 



This is the same equation as Eq. (|93[) . except that now 
the interesting region is < y < oo, and the "initial 
condition" is that \I/(oo) = 0. We define the integral 



1= / Vdy, 
lo 



(110) 
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FIG. 8: A comparison of methods for handling the super-Lya 
n > 3 two-photon decays. We show the numerical radiative 
transfer method (solid line) and the analytic approximation 
of Eqs. (flOO)) and (flTTj) (dashed line). 

which is a function of W. Then Eq. (j!03[) can be re- 
written as 



x 2 7 _ S ™Ly a T r ( X 2p ^h^ ya /T r 

+ c 3 n H h \3x ls 



(111) 



The integral X function of W is plotted in Fig. [6] 
We evaluate it with an implicit ODE each time it is 
needed. From Eq. (| 1 09[) we see that one should have 
1 -> as W -> 0, and * -> e~ v , 1 -> 1 as W -> oo. 
Figure [6] shows that this is precisely what happens. 

Equation (|100[) presents an implementation difficulty: 
it depends on the derivative of x 2 ^ , which itself depends 
on x e because of the implicit dependence of a?i s and x 2p in 
Eq. (| 1 1 lj) on x e . In principle one could solve this by treat- 
ing Eq. (jlOOp as an implicit equation for xP . In prac- 
tice, however, the effect of the two-photon decays from 
highly excited states is only a small perturbation on the 
recombination problem, so we have calculated Eq. (jlll|) 
using the "standard" recombination history described in 
Sec. HU We then calculate ±+ 7 for use in Eq. (|100|) by nu- 
merically differentiating this precomputed table of x+ (t ) . 

The result is shown in Fig. [5] One can see that the 
analytic approximation described here not only captures 
the qualitative fact that recombination is first accelerated 
and then delayed, because x 2 + first increases and then 
decreases, but also captures the quantitative magnitude 
of the effect. The maximum of x 2 ^ is 0.0061 at z = 1360, 
and one can see from Fig. [5] that this is indeed when the 
two-photon decays with v > is-L ya switch from speeding 
up recombination to slowing it down. 

The abundance x 2 ? of nonthermal photons above Lya 
from two-photon decay using the analytic estimate is 
shown in Fig. [9] 



D. Raman scattering 

We now consider Raman scattering. Selection rules 
imply that Raman scattering to the Is level is only pos- 
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FIG. 9: The number of nonthermal photons above Lya per 
H atom. The long-dashed curve shows the analytic estimate 
for x 2 ? (Sec. IV C|l ; the short-dashed curve shows the analytic 
estimate for x+ (Sec. |VT3)) ; and the solid curve shows their 
sum X 2 ? + a;5 . 



sible from s and d levels. The most populated such level 
is 2s, so we expect that Raman scattering is dominated 
by the 2s <-» Is process. This process contains contribu- 
tions from two types of photons: the photons near the 
peak of the blackbody spectrum, and the photons near 
Ha which contribute resonant Raman scattering, 



H(2s) + hv 1 -> H(3p) -► H(ls) + hu, 



(112) 



which is resonant when u' f=a and v m Vhyp- These 
pieces are (mostly) distinct: the radiation temperature 
is typically 0.02h7Z (z = 1200) whereas the Ha energy is 
Q.lAhTZ. Thus the rate of Raman scattering as a function 
of frequency v' of the incoming photon has two peaks - a 
nonresonant peak due to the large abundance of photons 
with energies of a few times T r , and a resonant peak at 
Ha. There are additional peaks due to the other Balmer 
resonances H/3, Hy, etc., but the abundance of these pho- 
tons is at least a factor of a few less than Ha (due to the 
Boltzmann factor) and the dipole matrix elements are 
smaller. This behavior can be seen in Fig. [T] 

The two sub-classes of Raman scattering transitions - 
nonresonant and resonant (Ha) - are considered sepa- 
rately here. 



1. Nonresonant scattering 

In terms of effects on recombination, the nonresonant 
Raman scattering is similar to the super- Lya two-photon 
decays in that a Raman scattering that produces an H I 
atom in the Is state also produces a photon above the 
Lya frequency. This photon will eventually redshift into 
the Lya line and excite a hydrogen atom, so again the net 
number of decays to the H I Is state is zero. Again the 
importance for recombination derives from the delay be- 
tween emission of the photon and its re-absorption. Thus 



Eq. (|100|) applies to Raman scattering as well, except that 
in the case of Raman scattering we will need a formula 



for 



the number of non-thermal photons above Lya 



contributed by nonresonant Raman scattering: 



.(Dl) 



(113) 



Our challenge now is to find an equation for x5. Re- 
viewing the steps that lead to Eq. (| 1 08[) shows that these 
are valid in the case of Raman scattering, except that 
we must replace the two-photon opacity with the Raman 
opacity, 



c 3 n H 
8nHv 3 



E 



9nl / 1 j. \ 
(1 + U'jXU 

9is 



fv'X2s 



dK n , 
dv 



(114) 

in place of Eq. |83|) . Taking the 2s term to be domi- 
nant, recalling that xi s <C x\ s , and using the blackbody 
formula for /„< , we can re- write this as 



c 3 n H 



dK? 



8ttHv 3 1 - e~y dv 



(115) 



where y = hv' /T^ = h(v — ti,y Q )/T r . 

The Raman scattering rate coefficient dK2 S /dv can be 
re- written using Eq. (|B10[) by noting that so long as v' <C 
vn ai we may take 2p to be the dominant intermediate 
state. Using the matrix element (2s||r||2p) 2 = 27a 2 ,: 

dK 2s o&:* 



dv 



\{U\\v\\2p)\ 



(116) 



The squared matrix element is related to the Lya decay 
rate via the dipole emission formula 



V 3 

fs TZ 2 a 2 



|(ls||r||2p)| 



from which we may derive 
dK 2s 9 



dv 



7T- ^Lya ( 



(117) 



(118) 



This is just the lowest-order term in the power series 
expansion of dKi s j dv\ there are corrections of order j/ 2 , 
etc. Due to the large prcfactor of the next-lowest 



,'•'5 



term we include it 
dK 2s 



dv 



— At (JL. 



-n 2 



1 + 8.15 



(119) 



The coefficient 8.15 was determined by numerical differ- 
entiation of dK2 S /dv. It is positive because of construc- 
tive interference between the dominant 2s — > 2p — > Is 
Raman scattering pathway and the neighboring 2s — > 
3p — > Is pathway. 

Making this approximation and plugging into 
Eq. (fTT5|) . we find 



9c 3 rt H ^Lya Xls 3 

167r 2 Hvl ya l-e~y ais 



8.15 



V' 

W 



(120) 
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We may re-write this using the Lya optical depth, 
Eq. (H|, as 



a fs TL ya l — 



1 



The Raman scattering version of Eq. (|108[) is then 

y 



f =v(l + 8.15^y ; 



where 



7 = 



(121) 
(122) 

(123) 



The normalization V is given by Eq. (| 1 23[) . Combining 
everything gives 



C R = 27 C ( 3 ) (II) 3 Tim ( 1 + 22.0^- 
+ 2^ y ' \hcj riH V ^ 



-hu hyoc /T t 



(128) 



The same implementation difficulty arises with 
Eq. (|113p as arose with Eq. (jlOOp . namely that in order 

to compute Ai| D ^ one needs to know which in turn 
requires one to already know the recombination history. 
We solve this problem by taking x*t from the standard 
recombination history (Sec.|TT|, just as we did in Sec. IV Cl 



Equation (|122[) has the integral solution, 

*(v) = v 

where 



y'e-y 



1 



1 + 8.15— y'\ e- Tiy ' v '^h,' 



dy', 
(124) 



r(y,y') = V 



y 



1 - e~y 



1 



■ 15 M y " 



is the optical depth between y and y' . We may then write 
the dimensionless integral of the spectral distortion: 



J 



^(y)dy 
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dy' 



1 - e~y' 



t \ ry 



X i 1+8 V)I dye 



(126) 



If r were negligible then the inner integral would collapse 
to give y' and the outer integral would be dominated by 
the region y' ~ 2. In practice for V <C 1 (the maximum is 
0.046 at z = 1184) r will be 1 in this region; for much 
larger y' it may be significant but all it does is suppress 
a portion of the integral that is already negligible. Thus 

we may neglect r and replace the integral dy e~ T ^ v ' v > 
with y' . This gives 



J 



V I dy' 
lo 



2C(3) + 8.15 
2((3)V [1 + 22.0 



y' 2 e~ v ' 
1 - e~y' 

T r tt 4 " 1 
MZ15 
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(127) 



We may now assemble all of the pieces and plug them 
into Eq. (|102[) : the spectral distortion is given in terms 
of * by Eq. fTDTj) . and the integral of * is Eq. (fT2T)) . 



2. Ha resonance 

Resonant Raman scattering via the 3p intermediate 
state [Eq. (|112p ] will induce a red damping wing correc- 
tion to the Ly/3 escape probability similar to the correc- 
tions discussed in Sec. lVBl for the two-photon corrections 
to Lya. Precisely the same calculation applies here ex- 
cept that the line width, Eq. (j9"Tj) gets replaced by 



r Ly/3 _ 

4tt 2 1 



.4 



3p.2s 



-hU3p,2s/T T 



(129) 



[The only differences in the calculation are the replace- 
ment of the initial state, 2s instead of nl; the intermediate 
state, 3p instead of 2p; the use of the Raman scattering 
rate coefficient dK/dv in place of dA/dv; since the low- 
frequency photon is in the initial state the use of l + f u ' in 
Eq. (|114|) instead of /„/ which appears in Eq. (|83"1) ; and 
no factor of (21 + l)/3 because the Einstein coefficient 
^3 P ,2s in Eq. Q129p has the states reversed from A n i t 2 P in 
Eq. (|9T|) .] Once again we can construct a dimensionless 
variable Wp = hWp/% and introduce a correction 



Az{ D2) 



A 



3p,ls 
TLy/3 



X 3 p[3>(-00) - 1], 



(130) 



where $(— oo) is evaluated from Eq. (f9"3")l using Wp in- 
stead of W. 

The parameter Wp measures the normalization of the 
Raman scattering opacity [to the Ly/3 photon, not the 
Ha photon; i.e. for the reverse reaction of Eq. (|112|) ] . 
just as W measured the normalization of the two-photon 
opacity to photons near Lya. We find that Wp reaches 
a maximum of 0.94 at z = 981. It is much larger than 
the corresponding W for the Lya line because two-photon 
absorption of a photon near Lya requires a three-particle 
interaction between the Lya photon, an atom, and a rare 
Balmer photon in the Wien tail of the CMB; whereas a 
photon near Ly/3 can be Raman-scattered in a two-body 
interaction without the help of any initial-state Balmer 
photons. The larger Wp means that the corrections to 
Ly/3 are of similar overall importance to the corrections 



22 



D. Effect from Raman scattering 



Error in analytic treatment of two-photon decays 
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FIG. 10: A comparison of methods for handling Raman scat- 
tering. We show the numerical radiative transfer method 
(solid line) and the analytic approximation of Eq. (I131|l and 
(fT28l) (dashed line). 



to Lya, even though in an absolute sense Ly/3 is less 
important. 

The overall change in the downward decay rate from 
Raman scattering is the sum of the nonrcsonant and res- 
onant contributions: 



Ax 



(Dl) 
I 



Ax 



(D2) 
i 



(131) 



However it should be noted that the downward rate D2 
(Eq. 1 1 30|) is emitted in the red damping wing of Ly/3 and 
is added to the Ly/3 spectral distortion fhyp- = f3p,is- 
appearing on the left hand side of Eq. ([6]). It is also 
included as an increase in the 3p — > Is transition rate 
x 3p ->i s . It therefore feeds back at later times on Lya. 
In contrast the rate Dl has all radiative transfer effects 
already included, and it is included along with Eq. (|100|) 
in the 2p — > Is rate. (It is not obvious whether it is better 
to analytically approximate Dl as being a net rate from 
the 2s or 2p state, since the Raman scattering occurs 
from 2s but re-excitations are to 2p. In practice 2s and 
2p are in equilibrium in the redshift range considered so 
this choice does not matter.) 

The effect of this correction on recombination is shown 
in Fig. 1101 Our analytic approximation reproduces most 
of the features of the full numerical result. The main 
deficiency is that the effects on Ax e /x e occur earlier in 
the analytic calculation than in the numerical result by 
about Az ~ 50. This error might be due to the approx- 
imation that x+ instantaneously reaches its steady-state 
value, but we defer investigation of this possibility to fu- 
ture work. 

As a summary, we show in Fig.[Tl]the total error in us- 
ing the analytical approximations presented here in place 
of more accurate numerical two-photon radiative trans- 
fer. The analytic approximations have reduced the max- 
imal fractional error in x e from 1.3% with no treatment 
of the two-photon effects to 0.3%, the maximal difference 
between the analytic and numerical treatments. 
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FIG. 11: The error in using analytic approximations for the 
two-photon effects. The vertical axis plots Ax e /x e , with the 
sign convention positive where the analytic approximation 
overestimates x e and negative where it underestimates x e . 



VI. EFFECT ON THE CMB 

We have calculated the effect of our corrected re- 
combination history on the CMB anisotropics using the 
CMBFast [45[ Boltzmann code. The results are shown 
in Fig. [T3J For the temperature and i?-mode polar- 
ization, we show the fractional change in the power 
spectrum Cj T or Cf E . For the cross-spectrum we 
show the change in correlation coefficient, defined by 

pTE = C TE^ C TT C EEy/2^ Thig ig mQre usefm than 

ACf E /Cf E , which is ill-behaved when Cj E passes 
through zero. For comparison we show the analytic ap- 
proximations of Sec. [V] (Note that we denote the CMB 
multipoles by £ to avoid confusion with atomic orbital 
angular momentum I.) 

While the mapping from the recombination history 
to the CMB power spectrum is in general complicated 
[§, [ll|, the qualitative features of Fig. [T^are easily un- 
derstood. The major effect on the temperature is that 
the decrease in free electron abundance in the 1100 < 
z < 1500 range increases the mean free path of the pho- 
tons and hence the Silk damping length [46l.l47l.l48l.l49t. 
This means that the oscillations of the baryon-photon 
fluid are damped more effectively when two-photon tran- 
sitions are included than in the standard picture, and 
hence the power spectrum is decreased. This effect is 
most significant at the higher multipoles where damping 
is most effective. In principle, changes in the recombi- 
nation history could also have moved the surface of last 
scatter, thereby shifting the acoustic peak positions and 
leading to oscillations in ACj T / Cf T ; however as one can 
see from the figure, this is a small effect. This is because 
the accelerating effect on recombination of sub-Lya two- 
photon decays and the decelerating effect of super-Lya 
decays and Raman scattering cancel at z k 1100. Thus, 
purely by accident, the surface of last scatter is not sig- 
nificantly shifted. 

Silk damping also suppresses the -E-mode polarization, 
and again this suppression is evident in the figure at 
high multipoles. However this is not the only effect of 
the new recombination history. In our calculation, re- 
combination is proceeding more slowly at the surface of 
last scatter, i.e. x e is decreasing more slowly than in the 
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standard picture, so the surface of last scatter is broader. 
Since generation of polarization requires a sufficiently low 
Thomson opacity to generate a local photon quadrupole 
by free streaming, and sufficiently high Thomson opacity 
to convert this quadrupole into polarization, the broader 
surface of last scatter enhances the -E-mode polarization 
p!^ ] . Since Silk damping increases rapidly as one goes to 
small scales, the broadening of the surface of last scatter 
wins at low I (£ < 1000) but Silk damping wins at high 
£, leading to the characteristic "positive, then negative" 
shape seen in the middle panel of Fig. [T2] At very low 
multipoles £ < 50 the old and new calculations agree de- 
spite the broader surface of last scatter; this is because 
the large-scale polarization is dominated by scattering 
at reionization, which we have left unchanged. However, 
the small changes described here are unobservable at such 
low £ because of cosmic variance. 

The temperature-polarization cross-correlation coeffi- 
cient is only slightly affected (< 0.002) because to a 
first approximation the Silk damping and the broadening 
of the last scattering surface only give a wavenumbcr- 
dependent re-scaling of the T and E. The only recom- 
bination process that would substantially change pJ E 
would be to move the surface of last scatter, and as noted 
above our corrections to recombination do not signifi- 
cantly move the surface of last scatter. 

One can estimate the importance of the changes in the 
recombination history for various CMB experiments by 
computing 



Change in power spectra 



(132) 



where Fui is the Fisher matrix for the temperature power 
spectra. If all cosmological parameters were known per- 
fectly and there were no systematic errors then Z repre- 
sents the number of sigma at which the changes in the 
recombination history could be detected. In the realistic 
case in which one is trying to measure the cosmologi- 
cal parameters, then the maximum possible bias from 
using the incorrect recombination history is Z sigma; 
the actual bias will depend on whether the parameter 
in question produces a similar or orthogonal change in 
the Ce's to the change from the recombination history. 
For the WMAP 5-year Fisher matrix F u < [3, H3] we find 
Z = 0.40, i.e. a 0.40cr effect for statistical errors only. 
In practice the calibration, beam, and point source un- 
certainties are significant in the range 200 < £ < 800 
relevant for WMAP [HI, [H, [H, Hi] ; if one includes them 
by adding these modes to the inverse Fisher matrix, this 
is reduced to 0.27c. A similar calculation is possible for 
the ACBAR final data release (55[; here one finds a 0.48a 
effect for the Fisher matrix only, which is reduced to 
0.19cr after marginalizing over the beam and calibration 
uncertainty. Thus the corrections described here have 
negligible impact on the interpretation of the currently 
published WMAP and ACBAR data @, H. However 
for the upcoming Planck mission [56| , if we use the data 
model from the Dark Energy Task Force (57| including 
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FIG. 12: The change in CMB temperature power spectrum, 
_E-mode polarization, and cross correlation due to the two- 
photon processes described in this paper. The full calculation 
is shown with the solid lines, and the prediction from our 
simple analytic model is shown with the dashed lines. 



temperature anisotropics in the 100, 143, and 217 GHz 
bands, and E-modc polarization in the 143 and 217 GHz 
bands [including Cf E in the sum and Fisher matrix in 
Eq. we find Z = 7, i.e. a 7cr effect. 

We have also calculated the importance of the error 
using the analytic approximation, i.e. we obtain Z 2 from 
Eq. (|132|) but using the difference between numerical and 
analytic results for ACj T . This gives Z = 0.8 for Planck, 
implying that the analytic treatment of two-photon de- 
cays in this paper is sufficient to reduce the errors from 7rj 
to 0.8(7 for Planck. While this is an enormous improve- 
ment, one would like to do better in the future, perhaps 
by computing the finite lag time required for x+ to reach 
its steady-state value, or (less attractively) with a fudge 
factor fit to the full numerical result. 

In addition to the CMB anisotropics, one could also 
consider the global spectral distortion from recombina- 
tion. While observation of this signal will be challeng- 
ing due to foregrounds [l6| and the sensitivity require- 
ments [HI , it has prompted a large amount of work over 
the years because it is a direct pro be of recombination 
physics [H, El, El, HI, HO, l6ll. l62l]. Our analysis has 
shown a ~ 1 order of magnitude increase in the spectral 



24 



distortion blucward of Lya; however all of these photons 
are absorbed and re-processed when they redshift into 
the Lya line and in our calculation the effect on the ob- 
servable spectral distortion (red side of Lya) is at the 
level of a few percent. We have not followed the lower- 
frequency distortion due to Balmer, Paschen, etc. lines 
in this paper, but the excited level occupations typically 
change by a few percent or less and we expect the final 
spectral distortion due to these lines to change by a sim- 
ilar amount. Given that the spectral distortion has not 
yet been detected, it seems unlikely that the corrections 
considered in this paper will be important for spectral 
distortion studies in the forseeable future. 



VII. CONCLUSION 

In this paper we have reconsidered the physics of two- 
photon decays in hydrogen recombination, with particu- 
lar emphasis on decays from highly excited levels n > 3 
of H i. We reviewed the conceptual difficulties that have 
plagued earlier attempts to incorporate these transitions, 
particularly those that arise from sequential "1+1" de- 
cays such as 3d — > 2p — > Is. We argued that the notion 
of an effective two-photon decay coefficient, useful and 
beloved for the 2s — ► Is decay of H I, must be aban- 
doned for the higher states in order to resolve these con- 
ceptual difficulties. Rather one must resolve these two- 
photon decays as a function of frequency of the emitted 
photons and solve the full radiative transfer problem for 
two-photon emission and absorption. 

In addition to the two-photon decays from the highly 
excited states, we have also considered Raman scatter- 
ing and two-photon recombination. These are physically 
and mathematically very similar processes, with the dif- 
ferences being moving one photon from the final to ini- 
tial state (for Raman scattering) and moving the initial 
electron from a bound to a free state (for two-photon 
recombination). Neither process has previously been in- 
vestigated in the cosmological context. 

Our major phcnomcnological findings are as follows: 

1. The traditional treatment of two-photon decays 
from 2s does not account for stimulated emission 
and re-absorption of the Lya spectral distortion. 
The maximal correction to x e from these effects is 
0.6% and over most of the redshift range the net 
effect is to delay recombination. We agree with 
both the physical picture and the calculation of 
past studies of these effects [H, HH ■ 

2. Recombination is sped up by sub- Lya two-photon 
decays from higher excited states; the most im- 
portant initial state is 3d. Physically, inclusion of 
this effect involves a combination of two-photon 
emission and absorption which must be handled 
through a radiative transfer calculation; there is 
no relevant effective two-photon decay coefficient 



A^d analogous to that for 2s. We find a maximum 
change of 1.7% in x e . 

3. Super- Lya two-photon decays initially speed up 
recombination by supplying another route to the 
ground state of hydrogen. These photons produce 
a delayed feedback that slows down recombination 
at z < 1100 as they redshift down to Lya. The 
maximum effect on x e is 0.7%. The time depen- 
dence of the line profile plays a central role in this 
effect. 

4. CMB photons can Raman scatter from hydrogen 
atoms in the 2s level and insert the atom into the 
ground state. Again this initially speeds up recom- 
bination but produces photons blueward of Lya 
that produce delayed feedback. The maximum ef- 
fect on x e is 1.0%. 

5. Direct two-photon recombinations to the ground 
state of hydrogen are negligible. 

6. The net effect of the two-photon transitions on the 
CMB power spectrum is two-fold: (i) they suppress 
the temperature and polarization power spectra by 
up to 2% at t ~ 3000 due to increased Silk damping 
from the faster early stages of recombination; and 
(ii) they increase the degree-scale polarization by 
~0.7% due to the wider surface of last scatter. At 
"Fisher matrix" level the net effect is la for Planck. 

We have constructed analytic approximations to all 
of the important two-photon processes that reproduce 
qualitatively and quantitatively the numerical radiative 
transfer results. The error from using these analytic ap- 
proximations instead of doing full two-photon radiative 
transfer will be 0.8<r for Planck. This is reassuring: we 
anticipate that once Planck identifies the relevant por- 
tion of parameter space, this ~ 0.8er correction to x e {z) 
can be added as a fudge factor. However much work lies 
ahead in order to make a full analysis of Planck data 
computationally feasible: our code even with the an- 
alytic approximations still takes ~ 1 day to run on a 
desktop machine, due primarily to the small step size of 
A In a = 4.25 x 10 -5 . We chose this time step because 
of our method of following feedback: the excitation in 
the Is — np line requires us to look up the photon phase 
space density at a slightly earlier time on the red side of 
the Is — (n + l)p line, and the time step must be small 
enough that the interpolation table for f\ s ,(n+i)p— nas 
already been built. Iterative feedback methods, such as 
those used in Rcf. [2t| . can use much larger time steps 
and hence should be better suited to fast calculation. 

In the future, we plan to incorporate these devel- 
opments into a more complete hydrogen recombination 
model that would include effects not considered here. 
Some of these effects, such as increasing the number of 
shells n max and collisions, have previously received treat- 
ment in the literature [25j, but ultimately one must si- 
multaneously incorporate all important processes. An- 
other process that may interact with some of the effects 
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here is Lya diffusion and recoil, for which one would 
have to add Fokker-Planck terms to the radiative transfer 
equations described herein. There is some work on this 
problem [2^, [36|, l37l l38l and for the analogous problem in 
the He I 594 A line [26| ; however given the importance of 
two-photon transitions, their optically thick nature, and 
the > 1% correction from effects that owe their existence 
to the time dependence of the line profile, it is clear that 
we need to solve the full time-dependent radiative prob- 
lem including both Lya diffusion and two-photon opacity. 
This will be the subject of future work. A final goal will 
be to speed up the recombination calculation to produce 
a recombination code suitable for inclusion into Markov 
Chain Monte Carlo codes for cosmological parameter es- 
timation. 
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APPENDIX A: MATTER TEMPERATURE 

This appendix constructs a solution for the matter 
temperature evolution valid at early times. We can write 
a formal integral solution to Eq. (|15[) since the equation 
is (if we fix x e ) a first-order linear ODE in T m /T r . The 
usual procedure gives 



T 



j(i> 



where 



J = 



o 1 + J(t) 



e -m K '(t)dt, 



3(1 + /ho + x e )m e cH 
and K is the integral defined by 



K{t) = J H(t)[l + J(t)]dt, 



(Al) 



(A2) 



(A3) 



with K(t) = and K'(t) < 0. The initial value or "par- 
ticular solution" is unimportant since K(Q) — oo. Since 



f 



-Kit] 



K'(t)di = —1, we may rewrite Eq. (|A1|) as 



i 



o l + J(t) 



e - K ®K'(t)dt. 



(A4) 



This is still only a formal solution because we cannot in 
general solve for x e without knowing T m . However in the 
case J > 1 which occurs during recombination (but is 
violated at z < 500), then the exponential factor e^ Ki - t " > 



kills the integrand in Eq. 
we may substitute J(t) ft 



l4J unless t w t. In this case 
J(i), and for J > 1 we get 



T 



1 



1 



j(ty 



(A5) 



i.e. Eq. 



APPENDIX B: TWO-PHOTON TRANSITION 
RATES 

This appendix is concerned with calculating the tran- 
sition rates for spontaneous two-photon decay (nl — > Is), 
Raman scattering (nl — > Is), and spontaneous two- 
photon recombination to the ground state of hydrogen. 
Stimulated rates and reverse reactions (Is — * nl and 
two-photon ionization) are obtainable by insertion of the 
usual 1 + f v factor and by the principle of detailed bal- 
ance, respectively. Note that the two-photon selection 
rules allow these reactions only when the initial state has 
/ = or I = 2 (sec DG05). This also holds for two-photon 
recombination although in this case it is acceptable to 
sum the two cases because for the H + + e~ continuum 
the /-sub-levels are statistically populated. 

We will have to consider both discrete and continuous 
states here. The latter appear in summations over in- 
termediate states in two-photon decay and Raman scat- 
tering, and as initial states in two-photon recombina- 
tion. They are characterized by an imaginary princi- 
pal quantum number to agree with the usual formula, 
E e = —IZ/n 2 (we take the 3n > branch), or alterna- 
tively by their wave number 



h Tin 



(Bl) 



The conventional normalization of the unbound states 
is that the radial wavefunction R(r) [with rip(r, 8, <fi) — 
R( r )Yim(9, </>)] is chosen to oscillate between —1 and +1, 
i.e. we normalize to a sphere of radius X = 2. In this 
case, continuum states are formally separated in wave 
number by A A; = ir/X = ir/2. Therefore in all calcu- 
lations of summations over intermediate p-states states, 
one should in actual computation make the replacement 



all n 



dk, 



(B2) 



where we will write "all n" for emphasis, and n 
ik/(jj}/ 2 K) for continuum states. 



1. Two-photon decay 

The spontaneous two-photon decay rate for the process 
H(nZ) ->■ H(ls) + hv + hv' (B3) 



2G 



is given by 



dA 



atv 3 v' 3 



dv 108(2/ + 1)TZ ( 



■\M\' 



(B4) 



where af s ~ 1/137 is the fine structure constant, 1Z is the 
hydrogen Rydberg constant, and the matrix element is 



M 



V(n/||r||iVp>(7Vp||r||i S } 



all TV 



1 - N~ 2 - v/TZ 1 - N- 2 - v'/TZ 



(B5) 



(See Eq. 13 of Ref. [22|.) The summation over interme- 
diate states is carried out over both the discrete values 
of N and the continuous spectrum (N~ 2 < 0, restricted 
of course to p symmetry) . Note that Ai as defined here 
is dimensionlcss, and hence so is dk/dv. 
One can write the total two-photon rate: 



f (i- n -*)n dA 
k n i = — dv. 

J \l-n- 2 )n/2 dv 



(B6) 



The cutoff at half the maximum frequency avoids double 
counting, since one photon is always above this cut and 
the other below. 



2. Raman scattering 



Now we are interested in the reaction 
H(nZ) + hv' -> H(ls) + hv. 



(B7) 



This reaction is related to the two-photon decay rate by 
crossing symmetry. Therefore, whereas the two-photon 
decay rate was 



dVnl.ls 



(27) 



dv v " 108(2/ + l)ft 6 
the Raman scattering rate is 



\M\ 2 (l+ }„)(!+}„,), (B8) 



dT n , 



Is 



dv 



(Raman) 



a%v 3 v' 3 
108(2/+ rjw 



■\M\ 2 {1+U)U^ (B9) 



where the matrix element Ai is given by Eq. ()B5[) . It is 
convenient for us to define the coefficient 



dK 



nl 



0$V 3 V 13 



dv 108(2/ + l)ft 6 



\M\ S 



(B10) 



that contains the matrix element and kinematic factors 
but not the state of the radiation field. Then 



dT nl _i s dK ni 

— (Raman) = — — (1 + j v )j v <. 

dv dv 

3. Two-photon recombination 



(Bll) 



We also require the rate coefficient for two-photon re- 
combination. 



H(ls) + hv + hv', 



(B12) 



where the initial-state electron has kinetic energy E e , and 
the final-state photons satisfy hv + hv' = E e — Ei s . We 
define the coefficient a^(v, v')n e n p as the number of re- 
combinations per unit time per unit frequency (dv) in 
a gas of monoenergetic elecrons and stationary protons, 
and in the absence of a stimulating radiation field. The 
computation of this rate is similar to that of the two- 
photon decay, except that now the initial state is in the 
continuum (n 2 < 0). This coefficient is a modification of 
Eq. |TB3]) : 



■\M\ 2 , (B13) 



a W(v,v')= J2 



1=0,2 



108(2/ + 1)U ( 



where N e (E, I) is the mean number of electrons per nor- 
malization volume with angular momentum /, assuming 
unit electron density. This is given by 



N e (E, I) = 2tt(2/ + 1) J = tt(2/ + 1) 



h 2 x 

/uE e ' 



(B14) 



where X is the radius of the normalization sphere. 
(This is most easily proved by considering the decom- 
position of a plane wave with unit probability density, 
e ik r into spherical harmonics, integrating the Zth par- 
tial wave out to radius X, and taking the large-X limit 
Jq X r 2 \ji(kr)\ 2 dr — > X/2k 2 .) This normalization volume 
cancels the X -1 / 2 dependence of the matrix element M. 
from the initial-state wave function. For our choice of 
normalization, X = 2 and we have 



>y)= E 



1=0,2 



(B15) 



The matrix element Ai is given by Eq. (|B5[) . Once again 
we restrict v > v'. Note that the normalization of the 
free state forces ip(r) to have units of cm -1 instead of 
cm -3 ' 2 , so for two-photon recombination Ai has units of 
cm 1 / 2 . Thus cS 2 \v, v') has units of cm 3 . 



4. Computation of matrix elements 

The two-photon matrix elements Ai can in principle 
be computed by brute force from Eq. (|B5[) . This was the 
approach taken in CS08 for hydrogen and in Hirata & 
Switzer for helium. Due to the cancellations among 
the various intermediate states for large n, however, we 
have chosen instead to use the Green function method, 
which is expected to be numerically stable even for very 
large n. The Green function method is described in the 
recombination context by Eqs. (24-26) of Ref. [22j , and 
more generally by Refs. |63l. |64||. 
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